DAMASK_EICMD/processing/post/table2ang.py

#!/usr/bin/env python

import os,re,sys,math,string,damask,numpy
from optparse import OptionParser, Option

# -----------------------------
class extendableOption(Option):
# -----------------------------
# used for definition of new option parser action 'extend', which enables to take multiple option arguments
# taken from online tutorial http://docs.python.org/library/optparse.html
  
  ACTIONS = Option.ACTIONS + ("extend",)
  STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
  TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
  ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)

  def take_action(self, action, dest, opt, value, values, parser):
    if action == "extend":
      lvalue = value.split(",")
      values.ensure_value(dest, []).extend(lvalue)
    else:
      Option.take_action(self, action, dest, opt, value, values, parser)

def integerFactorization(i):
  
  j = int(math.floor(math.sqrt(float(i))))
  while (j>1 and int(i)%j != 0):
    j -= 1
  return j

def positiveRadians(angle):

  angle = math.radians(float(angle))
  while angle < 0.0:
    angle += 2.0*math.pi

  return angle


def getHeader(sizeX,sizeY,step):
  
  return [ \
  '# TEM_PIXperUM          1.000000', \
  '# x-star                0.509548', \
  '# y-star                0.795272', \
  '# z-star                0.611799', \
  '# WorkingDistance       18.000000', \
  '#', \
  '# Phase                 1', \
  '# MaterialName          Al', \
  '# Formula               Fe', \
  '# Info', \
  '# Symmetry              43', \
  '# LatticeConstants      2.870 2.870 2.870  90.000  90.000  90.000', \
  '# NumberFamilies        4', \
  '# hklFamilies           1  1  0 1 0.000000 1', \
  '# hklFamilies           2  0  0 1 0.000000 1', \
  '# hklFamilies           2  1  1 1 0.000000 1', \
  '# hklFamilies           3  1  0 1 0.000000 1', \
  '# Categories            0 0 0 0 0 ', \
  '#', \
  '# GRID: SquareGrid', \
  '# XSTEP: ' + str(step), \
  '# YSTEP: ' + str(step), \
  '# NCOLS_ODD: ' + str(sizeX), \
  '# NCOLS_EVEN: ' + str(sizeX), \
  '# NROWS: ' + str(sizeY), \
  '#', \
  '# OPERATOR: ODFsammpling', \
  '#', \
  '# SAMPLEID: ', \
  '#', \
  '# SCANID: ', \
  '#', \
  ]

# --------------------------------------------------------------------
#                                MAIN
# --------------------------------------------------------------------

parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """
Builds an ang file out of ASCII table.

""" + string.replace('$Id$','\n','\\n')
)


parser.add_option('--coords',           dest='coords', type='string', \
                                        help='label of coords in ASCII table')
parser.add_option('--eulerangles',      dest='eulerangles', type='string', \
                                        help='label of euler angles in ASCII table')
parser.add_option('--defgrad',          dest='defgrad', type='string', \
                                        help='label of deformation gradient in ASCII table')
parser.add_option('-n','--normal',      dest='normal', type='float', nargs=3, \
                                        help='normal of slices to visualize')
parser.add_option('-s','--size',        dest='size', type='float', nargs=3, \
                                        help='physical size of ang file')
parser.add_option('-u','--up',          dest='up', type='float', nargs=3,
                                        help='up direction of slices to visualize')
parser.add_option('-r','--resolution',  dest='res', type='float',
                                        help='scaling factor for resolution')
parser.add_option('--hexagonal',   dest='hex', action='store_true',
                                        help='use in plane hexagonal grid')
parser.add_option('-c','--center',      dest='center', type='float', nargs=3,
                                        help='center of ang file in cube, negative for center')
parser.set_defaults(coords = 'coords')
parser.set_defaults(eulerangles = 'eulerangles')
parser.set_defaults(defgrad = 'f')
parser.set_defaults(hexagonal = False)
parser.set_defaults(normal = [0.0,0.0,1.0])
parser.set_defaults(size = [1.0,1.0,0.0])
parser.set_defaults(up = [1.0,0.0,0.0])
parser.set_defaults(center = [-1.0,-1.0,-1.0])
parser.set_defaults(res = 1.0)
(options,filenames) = parser.parse_args()

datainfo = {
             'vector':     {'len':3,
                            'label':[]},
             'tensor':     {'len':9,
                            'label':[]}
           }
           
datainfo['vector']['label'].append(options.coords)
datainfo['vector']['label'].append(options.eulerangles)
datainfo['tensor']['label'].append(options.defgrad)

# ------------------------------------------ setup file handles ---------------------------------------  

files = []
if filenames == []:
  files.append({'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout})
else:
  for name in filenames:
    if os.path.exists(name):
      files.append({'name':name, 'input':open(name)})


# ------------------------------------------ loop over input files ---------------------------------------  

for file in files:
  if file['name'] != 'STDIN': print file['name']

  table = damask.ASCIItable(file['input'])                                  # open ASCII_table for reading
  table.head_read()                                                         # read ASCII header info

# --------------- figure out dimension and resolution 

  try:
    locationCol = table.labels.index('ip.x')                                # columns containing location data
    
  except ValueError:
    print 'no coordinate data found...'
    continue

  grid = [{},{},{}]
  while table.data_read():                                                  # read next data line of ASCII table
    for j in xrange(3):
      grid[j][str(table.data[locationCol+j])] = True                        # remember coordinate along x,y,z
  resolution = numpy.array([len(grid[0]),\
                            len(grid[1]),\
                            len(grid[2]),],'i')                             # resolution is number of distinct coordinates found
  dimension = resolution/numpy.maximum(numpy.ones(3,'d'),resolution-1.0)* \
              numpy.array([max(map(float,grid[0].keys()))-min(map(float,grid[0].keys())),\
                           max(map(float,grid[1].keys()))-min(map(float,grid[1].keys())),\
                           max(map(float,grid[2].keys()))-min(map(float,grid[2].keys())),\
                          ],'d')                                            # dimension from bounding box, corrected for cell-centeredness
  if resolution[2] == 1:
    dimension[2] = min(dimension[:2]/resolution[:2])

  N = resolution.prod()
  print '\t%s @ %s'%(dimension,resolution)
  

 # --------------- figure out columns to process
  active = {}
  column = {}
  values = {}

  head = []
  for datatype,info in datainfo.items():
    for label in info['label']:
      key = {True :'1_%s',
             False:'%s'   }[info['len']>1]%label
      if key not in table.labels:
        sys.stderr.write('column %s not found...\n'%key)
      else:
        if datatype not in active:     active[datatype] = []
        if datatype not in column:     column[datatype] = {}
        if datatype not in values:     values[datatype] = {}
        active[datatype].append(label)
        column[datatype][label] = table.labels.index(key)                   # remember columns of requested data
        values[datatype][label] = numpy.array([0.0 for i in xrange(N*datainfo[datatype]['len'])])
        
# ------------------------------------------ read value field ---------------------------------------  

  table.data_rewind()
  idx = 0
  while table.data_read():                                                  # read next data line of ASCII table
    for datatype,labels in active.items():                                  # loop over vector,tensor
      for label in labels:                                                  # loop over all requested curls
        begin = idx*datainfo[datatype]['len']
        end = begin + datainfo[datatype]['len']
        values[datatype][label][begin:end]= numpy.array(map(float,table.data[column[datatype][label]:
                                             column[datatype][label]+datainfo[datatype]['len']]),'d')
    idx+=1

  stepSize = 0.0
  for i in xrange(3): stepSize+=dimension[i]/resolution[i]/3.0/options.res
  print 'step size', stepSize
 
  if options.hexagonal: 
    stepSize0 = stepSize * math.sin(1.0/3.0*math.pi)
  else: 
    stepSize0 = stepSize

  print 'step Size in x direction', stepSize0

  angRes = int(options.size[0]/stepSize0),\
           int(options.size[1]/stepSize),\
           max(int(options.size[2]/stepSize),1)
  print 'resolution of ang file', angRes

  if options.hexagonal: 
    NpointsSlice = angRes[0]//2*(angRes[1]-1)+(angRes[0]-angRes[0]//2)*angRes[1]
  else:
    NpointsSlice = angRes[0]*angRes[1]

  z = numpy.array(options.normal,dtype='float')
  z = z/numpy.linalg.norm(z)
  x = numpy.array(options.up,dtype='float')
  x = x/numpy.linalg.norm(x)
  y = numpy.cross(z,x)
  x = numpy.cross(y,z)
  print 'x unit vector', x, 'with norm ', numpy.linalg.norm(x)
  print 'y unit vector', y, 'with norm ', numpy.linalg.norm(y)
  print 'z unit vector', z, 'with norm ', numpy.linalg.norm(z)
  Favg = damask.core.math.tensorAvg(values['tensor']['%s'%(options.defgrad)].\
                         reshape(resolution[0],resolution[1],resolution[2],3,3))

  coordTransform = numpy.array([x,y,z])
  print 'rotation matrix', coordTransform

  mySlice = numpy.zeros(NpointsSlice*3)
  eulerangles = values['vector']['%s'%options.eulerangles].reshape([3,N],order='F')

  offset = ((dimension - options.size)/2.0 + (dimension/angRes)/2.0)/options.res
  print 'offset', offset
 # offset =  numpy.array([0.5,0.5,0.5],dtype='float')/[float(options.res[0]),float(options.res[1]),float(options.res[2])]*[dimension[0],dimension[1],dimension[2]]
  for i in xrange(angRes[2]):
    idx = 0
    for j in xrange(angRes[0]):
      if options.hexagonal: 
        res1=angRes[1]-j%2
        #myOffset = offset +float(j%2)* numpy.array([0.0,0.5,0.0],dtype='float')/[float(options.res[0]),float(options.res[1]),float(options.res[2])]*[dimension[0],dimension[1],dimension[2]]
        myOffset = offset +float(j%2)* numpy.array([0.0,0.5*stepSize,0.0],dtype='float')
      else:
        res1=angRes[1]
        myOffset = offset
      for k in xrange(res1):
        mySlice[idx*3:idx*3+3] = numpy.dot(coordTransform,[j*stepSize0,k*stepSize,i*stepSize]+myOffset)
        #print mySlice[idx*3:idx*3+3]
        idx+=1
    mySlice = mySlice.reshape([3,NpointsSlice],order='F')
    indices=damask.core.math.math_nearestNeighborSearch(3,Favg,numpy.array(
      dimension,dtype='float'),NpointsSlice,N,mySlice,values['vector']['%s'%options.coords].reshape([3,N],order='F'))/27
    fileOut=open(os.path.join(os.path.dirname(name),os.path.splitext(os.path.basename(name))[0]+'_%s.ang'%(angRes[2]-i-1)),'w')
    for line in getHeader(angRes[0],angRes[1],angRes[2]):
      fileOut.write(line + '\n')
  
    # write data
    for idx in xrange(NpointsSlice):
      fileOut.write(''.join(['%10.5f'%positiveRadians(angle) for angle in eulerangles[:,indices[idx]]])+
           ' %10.5f %10.5f'%(mySlice[1,idx],mySlice[0,idx])+
           ' 100.0 1.0 0 1 1.0\n')

    fileOut.close()