223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

added debug messages. 

added background message for estimated time (+ grain info).
fixed functionality—finally works as intended (lump orientations within threshold).
This commit is contained in:
Philip Eisenlohr 2015-11-20 16:27:57 +00:00
parent 38c536b6b7
commit 0db4a7fad0
1 changed files with 107 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
processing/pre/geom_fromTable.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys,string,math import os,sys,string,math,types,time
import scipy.spatial, numpy as np import scipy.spatial, numpy as np
from optparse import OptionParser from optparse import OptionParser
import damask import damask
@ -80,12 +80,16 @@ parser.add_option('--crystallite',
dest = 'crystallite', dest = 'crystallite',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'crystallite index to be used [%default]') help = 'crystallite index to be used [%default]')
parser.add_option('--debug',
dest = 'debug', action = 'store_true',
help = 'output debug info')
parser.set_defaults(symmetry = [damask.Symmetry.lattices[-1]], parser.set_defaults(symmetry = [damask.Symmetry.lattices[-1]],
tolerance = 0.0, tolerance = 0.0,
degrees = False, degrees = False,
homogenization = 1, homogenization = 1,
crystallite = 1, crystallite = 1,
debug = False,
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -110,7 +114,8 @@ if options.axes != None and not set(options.axes).issubset(set(['x','+x','-x','y
(options.quaternion,4,'quaternion'), (options.quaternion,4,'quaternion'),
(options.microstructure,1,'microstructure'), (options.microstructure,1,'microstructure'),
][np.where(input)[0][0]] # select input label that was requested ][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale all angles to radians
threshold = np.cos(options.tolerance/2.*toRadians) # cosine of (half of) tolerance angle
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
@ -134,23 +139,27 @@ for name in filenames:
errors = [] errors = []
if not 3 >= coordDim >= 2: if not 3 >= coordDim >= 2:
errors.append('coordinates {} need to have two or three dimensions.'.format(options.coordinates)) errors.append('coordinates "{}" need to have two or three dimensions.'.format(options.coordinates))
if not np.all(table.label_dimension(label) == dim): if not np.all(table.label_dimension(label) == dim):
errors.append('input {} needs to have dimension {}.'.format(label,dim)) errors.append('input "{}" needs to have dimension {}.'.format(label,dim))
if options.phase != None and table.label_dimension(options.phase) != 1: if options.phase and table.label_dimension(options.phase) != 1:
errors.append('phase column {} is not scalar.'.format(options.phase)) errors.append('phase column "{}" is not scalar.'.format(options.phase))
if errors != []: if errors != []:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
continue continue
table.data_readArray([options.coordinates,label]+([] if options.phase == None else [options.phase])) table.data_readArray([options.coordinates] \
+ ([label] if isinstance(label, types.StringTypes) else label) \
+ ([options.phase] if options.phase else []))
if coordDim == 2: if coordDim == 2:
table.data = np.insert(table.data,2,np.zeros(len(table.data)),axis=1) # add zero z coordinate for two-dimensional input table.data = np.insert(table.data,2,np.zeros(len(table.data)),axis=1) # add zero z coordinate for two-dimensional input
if options.debug: damask.util.croak('extending to 3D...')
if options.phase == None: if options.phase == None:
table.data = np.column_stack((table.data,np.ones(len(table.data)))) # add single phase if no phase column given table.data = np.column_stack((table.data,np.ones(len(table.data)))) # add single phase if no phase column given
if options.debug: damask.util.croak('adding dummy phase info...')
# --------------- figure out size and grid --------------------------------------------------------- # --------------- figure out size and grid ---------------------------------------------------------
@ -185,58 +194,104 @@ for name in filenames:
microstructure = table.data[:,colOri] microstructure = table.data[:,colOri]
nGrains = len(np.unique(microstructure)) nGrains = len(np.unique(microstructure))
else: else:
colPhase = colOri + np.sum(dim) # column of phase data comes after orientation
index = np.lexsort((table.data[:,0],table.data[:,1],table.data[:,2])) # index of rank when sorting x fast, z slow
rank = np.argsort(index) # rank of index
KDTree = scipy.spatial.KDTree((table.data[:,:3]-mincorner) / delta) # build KDTree with dX = dY = dZ = 1 and origin 0,0,0
microstructure = np.zeros(N,dtype = 'uint32') # initialize empty microstructure # --- start background messaging
symQuats = [] # empty list of sym equiv orientations
bg = damask.util.backgroundMessage()
bg.start()
colPhase = -1 # column of phase data comes last
bg.set_message('sorting positions...')
index = np.lexsort((table.data[:,0],table.data[:,1],table.data[:,2])) # index of position when sorting x fast, z slow
bg.set_message('building KD tree...')
KDTree = scipy.spatial.KDTree((table.data[index,:3]-mincorner) / delta) # build KDTree with dX = dY = dZ = 1 and origin 0,0,0
statistics = {'global': 0, 'local': 0}
grain = -np.ones(N,dtype = 'int32') # initialize empty microstructure
orientations = [] # empty list of orientations
multiplicity = [] # empty list of orientation multiplicity (number of group members)
phases = [] # empty list of phase info phases = [] # empty list of phase info
nGrains = 0 # counter for detected grains nGrains = 0 # counter for detected grains
myRank = 0 # rank of current grid point existingGrains = np.arange(nGrains)
myPos = 0 # position (in list) of current grid point
tick = time.clock()
bg.set_message('assigning grain IDs...')
for z in xrange(grid[2]): for z in xrange(grid[2]):
for y in xrange(grid[1]): for y in xrange(grid[1]):
for x in xrange(grid[0]): for x in xrange(grid[0]):
if (myRank+1)%(N/100.) < 1: damask.util.croak('.',False) if (myPos+1)%(N/500.) < 1:
myData = table.data[index[myRank]] time_delta = (time.clock()-tick) * (N - myPos) / myPos
mySym = options.symmetry[min(int(myData[colPhase]),len(options.symmetry))-1] # select symmetry from option (take last specified option for all with higher index) bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'%(time_delta//3600,time_delta%3600//60,time_delta%60,myPos,N,nGrains))
myData = table.data[index[myPos]] # read data for current grid point
myPhase = int(myData[colPhase])
mySym = options.symmetry[min(myPhase,len(options.symmetry))-1] # select symmetry from option (take last specified option for all with higher index)
if inputtype == 'eulers': if inputtype == 'eulers':
o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians, o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians,
symmetry = mySym).reduced() symmetry = mySym)
elif inputtype == 'matrix': elif inputtype == 'matrix':
o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3).transpose(), o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3).transpose(),
symmetry = mySym).reduced() symmetry = mySym)
elif inputtype == 'frame': elif inputtype == 'frame':
o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3], o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3],
myData[colOri[1]:colOri[1]+3], myData[colOri[1]:colOri[1]+3],
myData[colOri[2]:colOri[2]+3], myData[colOri[2]:colOri[2]+3],
)).reshape(3,3), )).reshape(3,3),
symmetry = mySym).reduced() symmetry = mySym)
elif inputtype == 'quaternion': elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = myData[colOri:colOri+4], o = damask.Orientation(quaternion = myData[colOri:colOri+4],
symmetry = mySym).reduced() symmetry = mySym)
neighbors = KDTree.query_ball_point([x,y,z], 3) # search points within radius cos_disorientations = -np.ones(1,dtype='f') # largest possible disorientation
breaker = False closest_grain = -1 # invalid neighbor
for n in neighbors: # check each neighbor neighbors = np.array(KDTree.query_ball_point([x,y,z], 3)) # point indices within radius
if myRank <= rank[n] or table.data[n,colPhase] != myData[colPhase]: continue # skip myself, anyone further ahead (cannot yet have a grain ID), and other phases neighbors = neighbors[(neighbors < myPos) & \
for symQ in symQuats[microstructure[rank[n]]-1]: (table.data[index[neighbors],colPhase] == myPhase)] # filter neighbors: skip myself, anyone further ahead (cannot yet have a grain ID), and other phases
if (symQ*o.quaternion).asAngleAxis(degrees = options.degrees)[0] <= options.tolerance: # found existing orientation resembling me grains = np.unique(grain[neighbors]) # unique grain IDs among valid neighbors
microstructure[myRank] = microstructure[rank[n]]
breaker = True; break
if breaker: break
if microstructure[myRank] == 0: # no other orientation resembled me if len(grains) > 0: # check immediate neighborhood first
nGrains += 1 # make new grain ... cos_disorientations = np.array([o.disorientation(orientations[grainID],
microstructure[myRank] = nGrains # ... and assign to me SST = False)[0].quaternion.w \
symQuats.append(o.symmetry.equivalentQuaternions(o.quaternion.conjugated())) # store all symmetrically equivalent orientations for future comparison for grainID in grains]) # store disorientation per grainID
phases.append(myData[colPhase]) # store phase info for future reporting closest_grain = np.argmax(cos_disorientations) # find grain among grains that has closest orientation to myself
match = 'local'
myRank += 1 if cos_disorientations[closest_grain] < threshold: # orientation not close enough?
grains = existingGrains[np.atleast_1d( ( np.array(phases) == myPhase ) & \
( np.in1d(existingGrains,grains,invert=True) ) )] # check every other already identified grain (of my phase)
damask.util.croak('') if len(grains) > 0:
cos_disorientations = np.array([o.disorientation(orientations[grainID],
SST = False)[0].quaternion.w \
for grainID in grains]) # store disorientation per grainID
closest_grain = np.argmax(cos_disorientations) # find grain among grains that has closest orientation to myself
match = 'global'
if cos_disorientations[closest_grain] >= threshold: # orientation now close enough?
grainID = grains[closest_grain]
grain[myPos] = grainID # assign myself to that grain ...
orientations[grainID] = damask.Orientation.average([orientations[grainID],o],
[multiplicity[grainID],1]) # update average orientation of best matching grain
multiplicity[grainID] += 1
statistics[match] += 1
else:
grain[myPos] = nGrains # ... and assign to me
orientations.append(o) # store new orientation for future comparison
multiplicity.append(1) # having single occurrence so far
phases.append(myPhase) # store phase info for future reporting
nGrains += 1 # update counter ...
existingGrains = np.arange(nGrains)
myPos += 1
bg.stop()
bg.join()
if options.debug: damask.util.croak("{} seconds total.\n{} local and {} global matches.".format(time.clock()-tick,statistics['local'],statistics['global']))
# --- generate header ---------------------------------------------------------------------------- # --- generate header ----------------------------------------------------------------------------
@ -248,11 +303,11 @@ for name in filenames:
'homogenization': options.homogenization, 'homogenization': options.homogenization,
} }
damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))), damask.util.croak(['grid a b c: {}'.format(' x '.join(map(str,info['grid']))),
'size x y z: %s'%(' x '.join(map(str,info['size']))), 'size x y z: {}'.format(' x '.join(map(str,info['size']))),
'origin x y z: %s'%(' : '.join(map(str,info['origin']))), 'origin x y z: {}'.format(' : '.join(map(str,info['origin']))),
'homogenization: %i'%info['homogenization'], 'homogenization: {}'.format(info['homogenization']),
'microstructures: %i'%info['microstructures'], 'microstructures: {}'.format(info['microstructures']),
]) ])
# --- write header --------------------------------------------------------------------------------- # --- write header ---------------------------------------------------------------------------------
@ -270,28 +325,22 @@ for name in filenames:
] ]
config_header += ['<texture>'] config_header += ['<texture>']
for i,quats in enumerate(symQuats): for i,orientation in enumerate(orientations):
config_header += ['[Grain%s]'%(str(i+1).zfill(formatwidth)), config_header += ['[Grain%s]'%(str(i+1).zfill(formatwidth)),
'axes\t%s %s %s'%tuple(options.axes) if options.axes != None else '', 'axes\t%s %s %s'%tuple(options.axes) if options.axes != None else '',
'(gauss)\tphi1 %g\tPhi %g\tphi2 %g\tscatter 0.0\tfraction 1.0'%tuple(np.degrees(quats[0].asEulers())), '(gauss)\tphi1 %g\tPhi %g\tphi2 %g\tscatter 0.0\tfraction 1.0'%tuple(orientation.asEulers(degrees = True)),
] ]
table.labels_clear() table.labels_clear()
table.info_clear() table.info_clear()
table.info_append([ table.info_append([scriptID + ' ' + ' '.join(sys.argv[1:])])
scriptID + ' ' + ' '.join(sys.argv[1:]), table.head_putGeom(info)
"grid\ta {}\tb {}\tc {}".format(*info['grid']), table.info_append(config_header)
"size\tx {}\ty {}\tz {}".format(*info['size']),
"origin\tx {}\ty {}\tz {}".format(*info['origin']),
"homogenization\t{}".format(info['homogenization']),
"microstructures\t{}".format(info['microstructures']),
config_header,
])
table.head_write() table.head_write()
# --- write microstructure information ------------------------------------------------------------ # --- write microstructure information ------------------------------------------------------------
table.data = microstructure.reshape(info['grid'][1]*info['grid'][2],info['grid'][0]) table.data = grain.reshape(info['grid'][1]*info['grid'][2],info['grid'][0]) + 1 # offset from starting index 0 to 1
table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ') table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ')
#--- output finalization -------------------------------------------------------------------------- #--- output finalization --------------------------------------------------------------------------