279 lines
11 KiB
Python
Executable File
279 lines
11 KiB
Python
Executable File
#!/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()
|
|
|