diff --git a/processing/post/3Dvisualize.py b/processing/post/3Dvisualize.py index d728dbbb1..baf096fd0 100755 --- a/processing/post/3Dvisualize.py +++ b/processing/post/3Dvisualize.py @@ -344,6 +344,10 @@ parser.add_option('--scaling', dest='scaling', type='float', \ help='scaling of fluctuation [%default]') parser.add_option('-u', '--unitlength', dest='unitlength', type='float', \ help='set unit length for 2D model [%default]') +parser.add_option('--filenodalcoords', dest='filenodalcoords', type='string', \ + help='ASCII table containing nodal coords') +parser.add_option('--labelnodalcoords', dest='nodalcoords', type='string', \ + help='labels of nodal coords in ASCII table') parser.add_option('-l', '--linear', dest='linearreconstruction', action='store_true',\ help='use linear reconstruction of geometry [%default]') @@ -362,12 +366,15 @@ parser.set_defaults(scaling = 1.0) parser.set_defaults(undeformed = False) parser.set_defaults(unitlength = 0.0) parser.set_defaults(cell = True) +parser.set_defaults(filenodalcoords = '') +parser.set_defaults(labelnodalcoords = 'coord') parser.set_defaults(linearreconstruction = False) sep = {'n': '\n', 't': '\t', 's': ' '} (options, args) = parser.parse_args() - +if options.scaling !=1.0 and options.linearreconstruction: print 'cannot scale for linear reconstruction' +if options.scaling !=1.0 and options.filenodalcoords!='': print 'cannot scale when reading coordinate from file' options.separator = options.separator.lower() for filename in args: if not os.path.exists(filename): @@ -452,24 +459,35 @@ for filename in args: dim[2] = min(dim/res) else: dim[2] = options.unitlength + + if options.filenodalcoords: + mesh = numpy.zeros(((res[0]+1)*(res[1]+1)*(res[2]+1),3),'d') + mesh=mesh.reshape(res[0]+1,res[1]+1,res[2]+1,3) + filenodalcoords = open(options.filenodalcoords) + tablenodalcoords = damask.ASCIItable(filenodalcoords) + tablenodalcoords.head_read() + coord = tablenodalcoords.labels.index(options.labelnodalcoords+'.x') + i = 0 + while tablenodalcoords.data_read(): + mesh[i%(res[0]+1),(i//(res[0]+1))%(res[1]+1),(i//(res[0]+1)//(res[1]+1)) % (res[2]+1),:]=\ + [float(tablenodalcoords.data[coord]), + float(tablenodalcoords.data[coord+1]), + float(tablenodalcoords.data[coord+2])] + i += 1 - if options.undeformed: - defgrad_av = numpy.eye(3) else: - defgrad_av = damask.core.math.tensorAvg(numpy.reshape(values[:,column['tensor'][options.defgrad]: - column['tensor'][options.defgrad]+9], - (res[0],res[1],res[2],3,3))) - if options.linearreconstruction: - centroids = damask.core.mesh.deformed_linear(res,dim,defgrad_av, - numpy.reshape(values[:,column['tensor'][options.defgrad]: - column['tensor'][options.defgrad]+9], - (res[0],res[1],res[2],3,3))) - else: - centroids = damask.core.mesh.deformed_fft(res,dim,defgrad_av,options.scaling, - numpy.reshape(values[:,column['tensor'][options.defgrad]: - column['tensor'][options.defgrad]+9], - (res[0],res[1],res[2],3,3))) - ms = damask.core.mesh.mesh_regular_grid(res,dim,defgrad_av,centroids) + F = numpy.reshape(values[:,column['tensor'][options.defgrad]: column['tensor'][options.defgrad]+9],(res[0],res[1],res[2],3,3)) + if options.undeformed: + Favg = numpy.eye(3) + else: + Favg = damask.core.math.tensorAvg(F) + + if options.linearreconstruction: + centroids = damask.core.mesh.deformed_linear(res,dim,Favg,F) + else: + centroids = damask.core.mesh.deformed_fft(res,dim,Favg,options.scaling,F) + mesh = damask.core.mesh.mesh_regular_grid(res,dim,Favg,centroids) + fields = {\ 'tensor': {},\ 'vector': {},\ @@ -505,9 +523,9 @@ for filename in args: print '\n' out = {} - if options.output_mesh: out['mesh'] = vtk_writeASCII_mesh(ms,fields,res,sep[options.separator]) + if options.output_mesh: out['mesh'] = vtk_writeASCII_mesh(mesh,fields,res,sep[options.separator]) if options.output_points: out['points'] = vtk_writeASCII_points(centroids,fields,res,sep[options.separator]) - if options.output_box: out['box'] = vtk_writeASCII_box(dim,defgrad_av) + if options.output_box: out['box'] = vtk_writeASCII_box(dim,Favg) for what in out.keys(): print what diff --git a/processing/post/addDeformedConfiguration.py b/processing/post/addDeformedConfiguration.py index c2b0a0ae2..2b859f8bf 100755 --- a/processing/post/addDeformedConfiguration.py +++ b/processing/post/addDeformedConfiguration.py @@ -48,9 +48,12 @@ parser.add_option('-c','--coordinates', dest='coords', type='string',\ help='column heading for coordinates [%default]') parser.add_option('-d','--defgrad', dest='defgrad', type='string', \ help='heading of columns containing tensor field values') - +parser.add_option('-l', '--linear', dest='linearreconstruction', action='store_true',\ + help='use linear reconstruction of geometry [%default]') + parser.set_defaults(coords = 'ip') parser.set_defaults(defgrad = 'f' ) +parser.set_defaults(linearreconstruction = False) (options,filenames) = parser.parse_args() @@ -115,7 +118,7 @@ for file in files: sys.stderr.write('column %s not found...\n'%key) else: defgrad = numpy.array([0.0 for i in xrange(N*9)]).reshape(list(res)+[3,3]) - table.labels_append(['%s_deformed'%(coord) for coord in 'x','y','z']) # extend ASCII header with new labels + table.labels_append(['ip_deformed.%s'%(coord) for coord in 'x','y','z']) # extend ASCII header with new labels column = table.labels.index(key) # ------------------------------------------ assemble header --------------------------------------- @@ -134,8 +137,10 @@ for file in files: # ------------------------------------------ process value field ---------------------------- defgrad_av = damask.core.math.tensorAvg(defgrad) - centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) - + if options.linearreconstruction: + centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) + else: + centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) # ------------------------------------------ process data --------------------------------------- table.data_rewind() diff --git a/processing/post/spectral_buildElements.py b/processing/post/nodesFromCentroids.py similarity index 71% rename from processing/post/spectral_buildElements.py rename to processing/post/nodesFromCentroids.py index b4959a637..4a761e7fd 100755 --- a/processing/post/spectral_buildElements.py +++ b/processing/post/nodesFromCentroids.py @@ -38,18 +38,21 @@ def index(location,res): # -------------------------------------------------------------------- parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """ -Calculates current coordinates and nodal displacement from IP/FP based deformation gradient. +Calculates current nodal coordinates and nodal displacement from IP/FP based data. """ + string.replace('$Id$','\n','\\n') ) -parser.add_option('-c','--coordinates', dest='coords', type='string',\ - help='column heading for coordinates [%default]') -parser.add_option('-d','--defgrad', dest='defgrad', type='string', \ - help='heading of columns containing tensor field values') +parser.add_option('-d','--deformed', dest='deformedCentroids', type='string',\ + help='column heading for position of deformed center of element/IP [%default]') +parser.add_option('-u','--undeformed', dest='undeformedCentroids', type='string',\ + help='column heading for position of undeformed center of element/IP [%default]') +parser.add_option('-f','--defgrad', dest='defgrad', type='string',\ + help='column heading for position of deformation gradient [%default]') -parser.set_defaults(coords = 'ip') -parser.set_defaults(defgrad = 'f' ) +parser.set_defaults(deformedCentroids = 'ip_deformed') +parser.set_defaults(undeformedCentroids = 'ip') +parser.set_defaults(defgrad = 'f') (options,filenames) = parser.parse_args() @@ -75,7 +78,7 @@ for file in files: # --------------- figure out dimension and resolution try: - locationCol = table.labels.index('%s.x'%options.coords) # columns containing location data + locationCol = table.labels.index('%s.x'%options.undeformedCentroids) # columns containing location data except ValueError: print 'no coordinate data found...' continue @@ -108,37 +111,49 @@ for file in files: # --------------- figure out columns to process + key = '%s.x' %options.deformedCentroids + if key not in table.labels: + sys.stderr.write('column %s not found...\n'%key) + else: + columnDeformed = table.labels.index(key) + + key = '%s.x' %options.undeformedCentroids + if key not in table.labels: + sys.stderr.write('column %s not found...\n'%key) + else: + columnUndeformed = table.labels.index(key) + key = '1_%s' %options.defgrad if key not in table.labels: sys.stderr.write('column %s not found...\n'%key) else: - column = table.labels.index(key) - - + columnDefgrad = table.labels.index(key) # ------------------------------------------ read value field --------------------------------------- - defgrad = numpy.array([0.0 for i in xrange(N*9)]).reshape(list(res)+[3,3]) + deformedCentroids = numpy.array([0.0 for i in xrange(N*3)]).reshape(list(res)+[3]) + undeformedCentroids = numpy.array([0.0 for i in xrange(N*3)]).reshape(list(res)+[3]) + defgrad = numpy.array([0.0 for i in xrange(N*9)]).reshape(list(res)+[3,3]) table.data_rewind() table.data_read() - inc = table.data[table.labels.index('inc')] + inc = int(eval(table.data[table.labels.index('inc')])) table.data_rewind() idx = 0 while table.data_read(): # read next data line of ASCII table (x,y,z) = location(idx,res) # figure out (x,y,z) position from line count idx += 1 - defgrad[x,y,z] = numpy.array(map(float,table.data[column:column+9]),'d').reshape(3,3) - + deformedCentroids[x,y,z] = numpy.array(map(float,table.data[columnDeformed:columnDeformed+3]),'d') + undeformedCentroids[x,y,z] = numpy.array(map(float,table.data[columnUndeformed:columnUndeformed+3]),'d') + defgrad[x,y,z] = numpy.array(map(float,table.data[columnDefgrad:columnDefgrad+9]),'d').reshape(3,3) file['input'].close() # close input ASCII table # ------------------------------------------ process value field ---------------------------- - - defgrad_av = damask.core.math.tensorAvg(defgrad) - centroids = damask.core.mesh.deformed_fft(res,geomdim,defgrad_av,1.0,defgrad) - nodes = damask.core.mesh.mesh_regular_grid(res,geomdim,defgrad_av,centroids) + defgrad_av = damask.core.math.math_tensorAvg(defgrad) + nodesDeformed = damask.core.mesh.mesh_regular_grid(res,geomdim,defgrad_av,deformedCentroids) + nodesUndeformed = damask.core.mesh.mesh_regular_grid(res,geomdim,numpy.identity(3,'d'),undeformedCentroids) # ------------------------------------------ process data --------------------------------------- @@ -148,7 +163,7 @@ for file in files: table.labels_append('inc elem node ip grain ') # extend ASCII header with new labels table.labels_append(['node.%s'%(coord) for coord in 'x','y','z']) # extend ASCII header with new labels - table.labels_append(['Displacement %s'%(coord) for coord in 'X','Y','Z']) # extend ASCII header with new labels + table.labels_append(['Displacement%s'%(coord) for coord in 'X','Y','Z']) # extend ASCII header with new labels table.head_write() ielem = 0 @@ -156,7 +171,7 @@ for file in files: for y in xrange(res[1]+1): for x in xrange(res[0]+1): ielem +=1 - entry = [inc,0,ielem,0,0,'\t'.join([str(a) for a in(nodes[x][y][z])]),'\t'.join([str(a) for a in (nodes[x][y][z] - (x,y,z)*(geomdim/res))])] + entry = [inc,0,ielem,0,0,'\t'.join([str(a) for a in(nodesUndeformed[x][y][z])]),'\t'.join([str(a) for a in (nodesUndeformed[x][y][z] - nodesDeformed[x][y][z])])] table.data_append(entry) table.data_write() # output processed line table.data_clear() diff --git a/processing/post/table2ang.py b/processing/post/table2ang.py new file mode 100755 index 000000000..1a058e9e2 --- /dev/null +++ b/processing/post/table2ang.py @@ -0,0 +1,194 @@ +#!/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.set_defaults(norm = 'ip') + + +(options,filenames) = parser.parse_args() + +datainfo = { + 'vector': {'len':3, + 'label':[]} + } + +datainfo['vector']['label'] += ['eulerangles'] +# ------------------------------------------ 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 = {} + divergence = {} + + 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] = {} + if datatype not in divergence: divergence[datatype] = {} + if label not in divergence[datatype]: divergence[datatype][label] = {} + 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() + print values['vector']['eulerangles'] + print numpy.shape(values['vector']['eulerangles']) + 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 + values[datatype][label][idx:idx+3]= numpy.array(map(float,table.data[column[datatype][label]: + column[datatype][label]+datainfo[datatype]['len']]),'d') + idx+=3 + for z in xrange(resolution[2]): + fileOut=open(os.path.join(os.path.dirname(name),os.path.splitext(os.path.basename(name))[0]+'_%s.ang'%z),'w') + for line in getHeader(resolution[0],resolution[1],1.0): + fileOut.write(line + '\n') + + # write data + for counter in xrange(resolution[0]*resolution[1]): + print z*resolution[0]*resolution[1]*3+counter*3,' - ',z*resolution[0]*resolution[1]*3+counter*3+3 + fileOut.write(''.join(['%10.5f'%positiveRadians(angle) for angle in values['vector']['eulerangles'][z*resolution[0]*resolution[1]*3+counter*3:z*resolution[0]*resolution[1]*3+counter*3+3]])+ + ''.join(['%10.5f'%coord for coord in [counter%resolution[0],counter//resolution[1]]])+ + ' 100.0 1.0 0 1 1.0\n') + counter += 1 + + fileOut.close() + diff --git a/processing/setup/setup_processing.py b/processing/setup/setup_processing.py index 91b216c8b..674c74249 100755 --- a/processing/setup/setup_processing.py +++ b/processing/setup/setup_processing.py @@ -116,8 +116,9 @@ bin_link = { \ 'postResults.py', 'spectral_iterationCount.py', 'spectral_parseLog.py', - 'spectral_buildElements.py', + 'nodesFromCentroids.py', 'tagLabel.py', + 'table2ang', ], }