From 7df8f04f65b33ef664abaab7f8066e1a57f78fe0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 4 Aug 2014 17:53:41 +0000 Subject: [PATCH] updated test for postprocessing and improve some of the scripts --- processing/post/addCauchy.py | 6 +- processing/post/addCompatibilityMismatch.py | 17 +- processing/post/addCurl.py | 4 +- processing/post/addDeformedConfiguration.py | 17 +- processing/post/addDeterminant.py | 31 ++-- processing/post/addDeviator.py | 37 ++-- processing/post/addDivergence.py | 4 +- processing/post/addEhkl.py | 36 ++-- processing/post/addEuclideanDistance.py | 2 +- processing/post/addIPFcolor.py | 27 +-- processing/post/addMises.py | 3 +- processing/post/addOrientations.py | 27 +-- processing/post/addPK2.py | 3 +- processing/post/addSpectralDecomposition.py | 149 +++++----------- processing/post/addStrainTensors.py | 188 ++++++++------------ 15 files changed, 226 insertions(+), 325 deletions(-) diff --git a/processing/post/addCauchy.py b/processing/post/addCauchy.py index a134b22ce..1c226e5a0 100755 --- a/processing/post/addCauchy.py +++ b/processing/post/addCauchy.py @@ -64,8 +64,7 @@ for file in files: for datatype,info in datainfo.items(): for label in info['label']: - key = {True :'1_%s', - False:'%s' }[info['len']>1]%label + key = '1_%s'%label if key not in table.labels: file['croak'].write('column %s not found...\n'%key) missingColumns = True # break if label not found @@ -77,7 +76,7 @@ for file in files: continue # ------------------------------------------ assemble header ------------------------------------ - table.labels_append(['%i_Cauchy'%(i+1) for i in xrange(datainfo['stress']['len'])]) # extend ASCII header with new labels + table.labels_append(['%i_Cauchy'%(i+1) for i in xrange(9)]) # extend ASCII header with new labels table.head_write() # ------------------------------------------ process data ---------------------------------------- @@ -87,7 +86,6 @@ for file in files: column['defgrad'][active['defgrad'][0]]+datainfo['defgrad']['len']]),'d').reshape(3,3) P = np.array(map(float,table.data[column['stress'][active['stress'][0]]: column['stress'][active['stress'][0]]+datainfo['stress']['len']]),'d').reshape(3,3) - table.data_append(list(1.0/np.linalg.det(F)*np.dot(P,F.T).reshape(9))) # [Cauchy] = (1/det(F)) * [P].[F_transpose] outputAlive = table.data_write() # output processed line diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index 1a9f1a8c7..83386a230 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -26,7 +26,7 @@ parser.add_option('--no-volume','-v', dest='noVolume', action='store_false', help='do not calculate volume mismatch [%default]') parser.add_option('-c','--coordinates', dest='coords', action='store', type='string', metavar='string', help='column heading for coordinates [%default]') -parser.add_option('-f','--deformation', dest='defgrad', action='store', type='string', metavar='string ', +parser.add_option('-f','--defgrad', dest='defgrad', action='store', type='string', metavar='string ', help='column heading for coordinates [%defgrad]') parser.set_defaults(noVolume = False) parser.set_defaults(noShape = False) @@ -82,14 +82,13 @@ for file in files: # --------------- figure out columns to process --------------------------------------------------- missingColumns = False - for datatype,info in datainfo.items(): - for label in info['label']: - key = '1_%s'%label - if key not in table.labels: - file['croak'].write('column %s not found...\n'%key) - missingColumns = True - else: - column = table.labels.index(key) # remember columns of requested data + for label in datainfo['defgrad']['label']: + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + missingColumns = True + else: + column = table.labels.index(key) # remember columns of requested data if missingColumns: continue diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index da33fd9bf..c7912918e 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -92,8 +92,7 @@ for file in files: for datatype,info in datainfo.items(): for label in info['label']: - key = {True :'1_%s', - False:'%s' }[info['len']>1]%label + key = '1_%s'%label if key not in table.labels: file['croak'].write('column %s not found...\n'%key) else: @@ -139,7 +138,6 @@ for file in files: for datatype,labels in active.items(): # loop over vector,tensor for label in labels: # loop over all requested norms table.data_append(list(curl[datatype][label][x,y,z].reshape(datainfo[datatype]['len']))) - outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- diff --git a/processing/post/addDeformedConfiguration.py b/processing/post/addDeformedConfiguration.py index fb340f7ba..0d0523255 100755 --- a/processing/post/addDeformedConfiguration.py +++ b/processing/post/addDeformedConfiguration.py @@ -23,7 +23,7 @@ Operates on periodic ordered three-dimensional data sets. parser.add_option('-c','--coordinates', dest='coords', action='store', type='string', metavar='string', help='column heading for coordinates [%default]') -parser.add_option('-d','--defgrad', dest='defgrad', action='store', type='string', metavar='string', +parser.add_option('-f','--defgrad', dest='defgrad', action='store', type='string', metavar='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]') @@ -80,14 +80,13 @@ for file in files: # --------------- figure out columns to process --------------------------------------------------- missingColumns = False - for datatype,info in datainfo.items(): - for label in info['label']: - key = '1_%s'%label - if key not in table.labels: - file['croak'].write('column %s not found...\n'%key) - missingColumns = True - else: - column = table.labels.index(key) # remember columns of requested data + for label in datainfo['defgrad']['label']: + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + missingColumns = True + else: + column = table.labels.index(key) # remember columns of requested data if missingColumns: continue diff --git a/processing/post/addDeterminant.py b/processing/post/addDeterminant.py index d9f8489e1..a6044c174 100755 --- a/processing/post/addDeterminant.py +++ b/processing/post/addDeterminant.py @@ -61,33 +61,28 @@ for file in files: table.head_read() # read ASCII header info table.info_append(string.replace(scriptID,'\n','\\n') + '\t' + ' '.join(sys.argv[1:])) - active = defaultdict(list) + active = [] column = defaultdict(dict) - 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: - file['croak'].write('column %s not found...\n'%key) - else: - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) # remember columns of requested data + for label in datainfo['tensor']['label']: + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + else: + active.append(label) + column[label] = table.labels.index(key) # remember columns of requested data # ------------------------------------------ assemble header --------------------------------------- - for datatype,labels in active.items(): # loop over vector,tensor - for label in labels: # loop over all requested determinants - table.labels_append('det(%s)'%label) # extend ASCII header with new labels + for label in active: + table.labels_append('det(%s)'%label) # extend ASCII header with new labels table.head_write() # ------------------------------------------ process data --------------------------------------- outputAlive = True while outputAlive and 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 determinantes - table.data_append(determinant(map(float,table.data[column[datatype][label]: - column[datatype][label]+datainfo[datatype]['len']]))) - + for label in active: + table.data_append(determinant(map(float,table.data[column[label]: + column[label]+datainfo['tensor']['len']]))) outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- diff --git a/processing/post/addDeviator.py b/processing/post/addDeviator.py index 72fb87d9d..b1b027817 100755 --- a/processing/post/addDeviator.py +++ b/processing/post/addDeviator.py @@ -65,36 +65,31 @@ for file in files: table.head_read() # read ASCII header info table.info_append(string.replace(scriptID,'\n','\\n') + '\t' + ' '.join(sys.argv[1:])) - active = defaultdict(list) + active = [] column = defaultdict(dict) - 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: - file['croak'].write('column %s not found...\n'%key) - else: - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) # remember columns of requested data + for label in datainfo['tensor']['label']: + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + else: + active.append(label) + column[label] = table.labels.index(key) # remember columns of requested data # ------------------------------------------ assemble header --------------------------------------- - for datatype,labels in active.items(): # loop over vector,tensor - for label in labels: # loop over all requested determinants - table.labels_append(['%i_dev(%s)'%(i+1,label) for i in xrange(9)]) # extend ASCII header with new labels - if(options.hydrostatic): table.labels_append('sph(%s)'%label) + for label in active: + table.labels_append(['%i_dev(%s)'%(i+1,label) for i in xrange(9)]) # extend ASCII header with new labels + if(options.hydrostatic): table.labels_append('sph(%s)'%label) table.head_write() # ------------------------------------------ process data --------------------------------------- outputAlive = True while outputAlive and 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 deviators - myTensor = map(float,table.data[column[datatype][label]: - column[datatype][label]+datainfo[datatype]['len']]) - table.data_append(deviator(myTensor)) - if(options.hydrostatic): table.data_append(oneThird*(myTensor[0]+myTensor[4]+myTensor[8])) - + for label in active: + myTensor = map(float,table.data[column[label]: + column[label]+datainfo['tensor']['len']]) + table.data_append(deviator(myTensor)) + if(options.hydrostatic): table.data_append(oneThird*(myTensor[0]+myTensor[4]+myTensor[8])) outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 950c92250..b2badb42b 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -106,8 +106,7 @@ for file in files: for datatype,info in datainfo.items(): for label in info['label']: - key = {True :'1_%s', - False:'%s' }[info['len']>1]%label + key = '1_%s'%label if key not in table.labels: file['croak'].write('column %s not found...\n'%key) else: @@ -161,7 +160,6 @@ for file in files: for label in labels: # loop over all requested for accuracy in options.accuracy: table.data_append(list(divergence[datatype][label][accuracy][x,y,z].reshape(datainfo[datatype]['len']//3))) - outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- diff --git a/processing/post/addEhkl.py b/processing/post/addEhkl.py index 6efcb2b9c..edbced28b 100755 --- a/processing/post/addEhkl.py +++ b/processing/post/addEhkl.py @@ -77,36 +77,30 @@ for file in files: table.head_read() # read ASCII header info table.info_append(string.replace(scriptID,'\n','\\n') + '\t' + ' '.join(sys.argv[1:])) - active = defaultdict(list) + active = [] column = defaultdict(dict) - for datatype,info in datainfo.items(): - for label in info['label']: - foundIt = False - for key in ['1_'+label,label]: - if key in table.labels: - foundIt = True - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) # remember columns of requested data - if not foundIt: - file['croak'].write('column %s not found...\n'%label) + for label in datainfo['vector']['label']: + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + else: + active.append(label) + column[label] = table.labels.index(key) # remember columns of requested data # ------------------------------------------ assemble header --------------------------------------- - for datatype,labels in active.items(): # loop over vector,tensor - for label in labels: # loop over all requested stiffnesses - table.labels_append('E%i%i%i(%s)'%(options.hkl[0], - options.hkl[1], - options.hkl[2],label)) # extend ASCII header with new labels + for label in active: + table.labels_append('E%i%i%i(%s)'%(options.hkl[0], + options.hkl[1], + options.hkl[2],label)) # extend ASCII header with new labels table.head_write() # ------------------------------------------ process data ---------------------------------------- outputAlive = True while outputAlive and 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 stiffnesses - table.data_append(E_hkl(map(float,table.data[column[datatype][label]:\ - column[datatype][label]+datainfo[datatype]['len']]),options.hkl)) - + for label in active: + table.data_append(E_hkl(map(float,table.data[column[label]:\ + column[label]+datainfo['vector']['len']]),options.hkl)) outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- diff --git a/processing/post/addEuclideanDistance.py b/processing/post/addEuclideanDistance.py index 3a75f6d80..f0acde8ff 100755 --- a/processing/post/addEuclideanDistance.py +++ b/processing/post/addEuclideanDistance.py @@ -194,8 +194,8 @@ for file in files: while table.data_read(): for i in xrange(len(feature_list)): table.data_append(distance[i,l]) # add all distance fields - outputAlive = table.data_write() # output processed line l += 1 + outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- outputAlive and table.output_flush() # just in case of buffered ASCII table diff --git a/processing/post/addIPFcolor.py b/processing/post/addIPFcolor.py index bf9d1de53..b5e0008d8 100755 --- a/processing/post/addIPFcolor.py +++ b/processing/post/addIPFcolor.py @@ -61,6 +61,11 @@ if options.a != None and \ if options.matrix != None: datainfo['tensor']['label'] += [options.matrix]; input = 'matrix' if options.quaternion != None: datainfo['quaternion']['label'] += [options.quaternion]; input = 'quaternion' +inputGiven = 0 +for datatype,info in datainfo.items(): + inputGiven += len(info['label']) +if inputGiven != 1: parser.error('select exactly one input format...') + toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians pole = np.array(options.pole) pole /= np.linalg.norm(pole) @@ -85,18 +90,20 @@ for file in files: active = defaultdict(list) column = defaultdict(dict) + missingColumns = False for datatype,info in datainfo.items(): for label in info['label']: - foundIt = False - for key in ['1_'+label,label]: - if key in table.labels: - foundIt = True - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) # remember columns of requested data - if not foundIt: - file['croak'].write('column %s not found...\n'%label) - break + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + missingColumns = True # break if label not found + else: + active[datatype].append(label) + column[datatype][label] = table.labels.index(key) # remember columns of requested data + + if missingColumns: + continue # ------------------------------------------ assemble header --------------------------------------- table.labels_append(['%i_IPF_%g%g%g'%(i+1,options.pole[0],options.pole[1],options.pole[2]) for i in xrange(3)]) @@ -124,7 +131,7 @@ for file in files: symmetry=options.symmetry).reduced() elif input == 'quaternion': o = damask.Orientation(quaternion=np.array(map(float,table.data[column['quaternion'][options.quaternion]:\ - column['quaternion'][options.quaternion]+datainfo['quaternion']['len']])), + column['quaternion'][options.quaternion]+datainfo['quaternion']['len']])), symmetry=options.symmetry).reduced() table.data_append(o.IPFcolor(pole)) diff --git a/processing/post/addMises.py b/processing/post/addMises.py index 8a58691b5..f6b5ca15b 100755 --- a/processing/post/addMises.py +++ b/processing/post/addMises.py @@ -75,8 +75,7 @@ for file in files: for datatype,info in datainfo.items(): for label in info['label']: - key = {True :'1_%s', - False:'%s' }[info['len']>1]%label + key = '1_%s'%label if key not in table.labels: file['croak'].write('column %s not found...\n'%key) else: diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index efb691007..c1e99f23c 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -69,7 +69,12 @@ if options.a != None and \ if options.matrix != None: datainfo['tensor']['label'] += [options.matrix]; input = 'matrix' if options.quaternion != None: datainfo['quaternion']['label'] += [options.quaternion]; input = 'quaternion' -toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians +inputGiven = 0 +for datatype,info in datainfo.items(): + inputGiven += len(info['label']) +if inputGiven != 1: parser.error('select exactly one input format...') + +toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians options.output = map(lambda x: x.lower(), options.output) r = damask.Quaternion().fromAngleAxis(toRadians*options.rotation[0],options.rotation[1:]) @@ -95,18 +100,20 @@ for file in files: active = defaultdict(list) column = defaultdict(dict) + missingColumns = False for datatype,info in datainfo.items(): for label in info['label']: - foundIt = False - for key in ['1_'+label,label]: - if key in table.labels: - foundIt = True - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) # remember columns of requested data - if not foundIt: - file['croak'].write('column %s not found...\n'%label) - break + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + missingColumns = True # break if label not found + else: + active[datatype].append(label) + column[datatype][label] = table.labels.index(key) # remember columns of requested data + + if missingColumns: + continue # ------------------------------------------ assemble header --------------------------------------- for output in options.output: diff --git a/processing/post/addPK2.py b/processing/post/addPK2.py index 8eb7ba0c2..3821af6ab 100755 --- a/processing/post/addPK2.py +++ b/processing/post/addPK2.py @@ -64,8 +64,7 @@ for file in files: for datatype,info in datainfo.items(): for label in info['label']: - key = {True :'1_%s', - False:'%s' }[info['len']>1]%label + key = '1_%s'%label if key not in table.labels: file['croak'].write('column %s not found...\n'%key) missingColumns = True # break if label not found diff --git a/processing/post/addSpectralDecomposition.py b/processing/post/addSpectralDecomposition.py index c373c66a4..2dbd10431 100755 --- a/processing/post/addSpectralDecomposition.py +++ b/processing/post/addSpectralDecomposition.py @@ -1,42 +1,27 @@ #!/usr/bin/env python # -*- coding: UTF-8 no BOM -*- -import os,re,sys,math,numpy,string,damask -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) +import os,re,sys,math,string +import numpy as np +from collections import defaultdict +from optparse import OptionParser +import damask +scriptID = '$Id$' +scriptName = scriptID.split()[1] # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- -parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """ +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ Add column(s) containing eigenvalues and eigenvectors of requested tensor column(s). -""" + string.replace('$Id$','\n','\\n') +""", version = string.replace(scriptID,'\n','\\n') ) - -parser.add_option('-t','--tensor', dest='tensor', action='extend', type='string', \ +parser.add_option('-t','--tensor', dest='tensor', action='extend', type='string', metavar='', help='heading of columns containing tensor field values') - parser.set_defaults(tensor = []) (options,filenames) = parser.parse_args() @@ -44,97 +29,59 @@ parser.set_defaults(tensor = []) if len(options.tensor) == 0: parser.error('no data column specified...') -datainfo = { # list of requested labels per datatype +datainfo = { # list of requested labels per datatype 'tensor': {'len':9, 'label':[]}, } +datainfo['tensor']['label'] += options.tensor -if options.tensor != None: datainfo['tensor']['label'] += options.tensor - -# ------------------------------------------ setup file handles --------------------------------------- - +# ------------------------------------------ setup file handles ------------------------------------ files = [] -if filenames == []: - files.append({'name':'STDIN', - 'input':sys.stdin, - 'output':sys.stdout, - 'croak':sys.stderr, - }) -else: - for name in filenames: - if os.path.exists(name): - files.append({'name':name, - 'input':open(name), - 'output':open(name+'_tmp','w'), - 'croak':sys.stdout, - }) - -# ------------------------------------------ loop over input files --------------------------------------- +for name in filenames: + if os.path.exists(name): + files.append({'name':name, 'input':open(name), 'output':open(name+'_tmp','w'), 'croak':sys.stderr}) +#--- loop over input files ------------------------------------------------------------------------ for file in files: - if file['name'] != 'STDIN': file['croak'].write(file['name']+'\n') + file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n') + table = damask.ASCIItable(file['input'],file['output'],True) # make unbuffered ASCII_table + table.head_read() # read ASCII header info + table.info_append(string.replace(scriptID,'\n','\\n') + '\t' + ' '.join(sys.argv[1:])) - table = damask.ASCIItable(file['input'],file['output'],False) # make unbuffered ASCII_table - table.head_read() # read ASCII header info - table.info_append(string.replace('$Id$','\n','\\n') + \ - '\t' + ' '.join(sys.argv[1:])) - - active = {} - column = {} - 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: - file['croak'].write('column %s not found...\n'%key) - else: - if datatype not in active: active[datatype] = [] - if datatype not in column: column[datatype] = {} - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) # remember columns of requested data - table.labels_append(['%i_eigval(%s)'%(i+1,label) - for i in xrange(3)]) # extend ASCII header with new labels - table.labels_append(['%i_eigvec(%s)'%(i+1,label) - for i in xrange(9)]) # extend ASCII header with new labels + active = [] + column = defaultdict(dict) + + for label in datainfo['tensor']['label']: + key = '1_%s'%label + if key not in table.labels: + file['croak'].write('column %s not found...\n'%key) + else: + active.append(label) + column[label] = table.labels.index(key) # remember columns of requested data # ------------------------------------------ assemble header --------------------------------------- - + for labels in active: + table.labels_append(['%i_eigval(%s)'%(i+1,label) for i in xrange(3)]) # extend ASCII header with new labels + table.labels_append(['%i_eigvec(%s)'%(i+1,label) for i in xrange(9)]) # extend ASCII header with new labels table.head_write() -# ------------------------------------------ process data --------------------------------------- - - 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 norms - tensor = numpy.array(map(float,table.data[column[datatype][label]: - column[datatype][label]+datainfo[datatype]['len']])).reshape((datainfo[datatype]['len']//3,3)) - (u,v) = numpy.linalg.eig(tensor) - table.data_append(list(u)) - table.data_append(list(v.transpose().reshape(datainfo[datatype]['len']))) - - table.data_write() # output processed line +# ------------------------------------------ process data ---------------------------------------- + outputAlive = True + while outputAlive and table.data_read(): # read next data line of ASCII table + for labels in active: # loop over requested data + tensor = np.array(map(float,table.data[column[label]:column[label]+datainfo['tensor']['len']])).\ + reshape((datainfo['tensor']['len']//3,3)) + (u,v) = np.linalg.eig(tensor) + table.data_append(list(u)) + table.data_append(list(v.transpose().reshape(datainfo['tensor']['len']))) + outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- + outputAlive and table.output_flush() # just in case of buffered ASCII table - table.output_flush() # just in case of buffered ASCII table - - try: - file['output'].close() # close output ASCII table - except: - pass - try: - file['croak'].close() # stop croaking - except: - pass - try: - file['input'].close() # close input ASCII table - except: - pass - + file['input'].close() # close input ASCII table (works for stdin) + file['output'].close() # close output ASCII table (works for stdout) if file['name'] != 'STDIN': - os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new + os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new diff --git a/processing/post/addStrainTensors.py b/processing/post/addStrainTensors.py index b57f38462..1554097fa 100755 --- a/processing/post/addStrainTensors.py +++ b/processing/post/addStrainTensors.py @@ -1,74 +1,54 @@ #!/usr/bin/env python # -*- coding: UTF-8 no BOM -*- -import os,re,sys,math,numpy,string,damask -from optparse import OptionParser, Option +import os,re,sys,math,string +import numpy as np +from collections import defaultdict +from optparse import OptionParser +import damask scriptID = '$Id$' scriptName = scriptID.split()[1] -# ----------------------------- -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 operator(stretch,strain,eigenvalues): - return { \ - 'V#ln': numpy.log(eigenvalues) , - 'U#ln': numpy.log(eigenvalues) , - 'V#Biot': ( numpy.ones(3,'d') - 1.0/eigenvalues ) , - 'U#Biot': ( eigenvalues - numpy.ones(3,'d') ) , - 'V#Green': ( numpy.ones(3,'d') - 1.0/eigenvalues*eigenvalues) *0.5, - 'U#Green': ( eigenvalues*eigenvalues - numpy.ones(3,'d')) *0.5, + return { + 'V#ln': np.log(eigenvalues) , + 'U#ln': np.log(eigenvalues) , + 'V#Biot': ( np.ones(3,'d') - 1.0/eigenvalues ) , + 'U#Biot': ( eigenvalues - np.ones(3,'d') ) , + 'V#Green': ( np.ones(3,'d') - 1.0/eigenvalues*eigenvalues) *0.5, + 'U#Green': ( eigenvalues*eigenvalues - np.ones(3,'d')) *0.5, }[stretch+'#'+strain] - # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- -parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """ +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ Add column(s) containing given strains based on given stretches of requested deformation gradient column(s). """ + string.replace(scriptID,'\n','\\n') ) - -parser.add_option('-u','--right', action='store_true', dest='right', \ - help='material strains based on right Cauchy--Green deformation, i.e., C and U') -parser.add_option('-v','--left', action='store_true', dest='left', \ - help='spatial strains based on left Cauchy--Green deformation, i.e., B and V') -parser.add_option('-l','-0','--logarithmic', action='store_true', dest='logarithmic', \ - help='calculate logarithmic strain tensor') -parser.add_option('-b','-1','--biot', action='store_true', dest='biot', \ - help='calculate biot strain tensor') -parser.add_option('-g','-2','--green', action='store_true', dest='green', \ - help='calculate green strain tensor') -parser.add_option('-f','--deformation', dest='defgrad', action='extend', type='string', \ - help='heading(s) of columns containing deformation tensor values [f]') - +parser.add_option('-u','--right', dest='right', action='store_true', + help='material strains based on right Cauchy--Green deformation, i.e., C and U [%default]') +parser.add_option('-v','--left', dest='left', action='store_true', + help='spatial strains based on left Cauchy--Green deformation, i.e., B and V [%default]') +parser.add_option('-l','-0','--logarithmic', dest='logarithmic', action='store_true', + help='calculate logarithmic strain tensor [%default]') +parser.add_option('-b','-1','--biot', dest='biot', action='store_true', + help='calculate biot strain tensor [%default]') +parser.add_option('-g','-2','--green', dest='green', action='store_true', + help='calculate green strain tensor [%default]') +parser.add_option('-f','--defgrad', dest='defgrad', action='extend', type='string', metavar = '', + help='heading(s) of columns containing deformation tensor values %default') parser.set_defaults(right = False) parser.set_defaults(left = False) parser.set_defaults(logarithmic = False) parser.set_defaults(biot = False) parser.set_defaults(green = False) -parser.set_defaults(defgrad = []) +parser.set_defaults(defgrad = ['f']) (options,filenames) = parser.parse_args() @@ -82,18 +62,14 @@ if options.logarithmic: strains.append('ln') if options.biot: strains.append('Biot') if options.green: strains.append('Green') -datainfo = { # list of requested labels per datatype +datainfo = { # list of requested labels per datatype 'defgrad': {'len':9, 'label':[]}, } -if options.defgrad == []: - datainfo['defgrad']['label'] = ['f'] -else: - datainfo['defgrad']['label'] = options.defgrad +datainfo['defgrad']['label'] = options.defgrad # ------------------------------------------ setup file handles --------------------------------------- - files = [] if filenames == []: files.append({'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout, 'croak':sys.stderr}) @@ -102,78 +78,68 @@ else: if os.path.exists(name): files.append({'name':name, 'input':open(name), 'output':open(name+'_tmp','w'), 'croak':sys.stderr}) -# ------------------------------------------ loop over input files --------------------------------------- - +# ------------------------------------------ loop over input files --------------------------------------- for file in files: if file['name'] != 'STDIN': file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n') else: file['croak'].write('\033[1m'+scriptName+'\033[0m\n') - table = damask.ASCIItable(file['input'],file['output'],False) # make unbuffered ASCII_table - table.head_read() # read ASCII header info + table = damask.ASCIItable(file['input'],file['output'],False) # make unbuffered ASCII_table + table.head_read() # read ASCII header info table.info_append(string.replace(scriptID,'\n','\\n') + '\t' + ' '.join(sys.argv[1:])) - active = {} - column = {} - head = [] + active = [] + column = defaultdict(dict) - 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] = {} - active[datatype].append(label) - column[datatype][label] = table.labels.index(key) - for theStretch in stretches: - for theStrain in strains: - table.labels_append(['%i_%s(%s)%s'%(i+1,theStrain,theStretch, - {True: label,False: ''}[label!='f']) - for i in xrange(datainfo['defgrad']['len'])]) # extend ASCII header with new labels + for label in datainfo['defgrad']['label']: + key = '1_%s'%label + if key not in table.labels: + sys.stderr.write('column %s not found...\n'%key) + else: + active.append(label) + column[label] = table.labels.index(key) -# ------------------------------------------ assemble header --------------------------------------- +# ------------------------------------------ assemble header --------------------------------------- + for label in active: + for theStretch in stretches: + for theStrain in strains: + table.labels_append(['%i_%s(%s)%s'%(i+1,theStrain,theStretch, + {True: label,False: ''}[label!='f'])for i in xrange(9)]) # extend ASCII header with new labels table.head_write() -# ------------------------------------------ process data --------------------------------------- +# ------------------------------------------ process data ---------------------------------------- + outputAlive = True + while outputAlive and table.data_read(): # read next data line of ASCII table + for label in active: # loop over all requested norms + F = np.array(map(float,table.data[column[label]: + column[label]+datainfo['defgrad']['len']]),'d').reshape(3,3) + (U,S,Vh) = np.linalg.svd(F) + R = np.dot(U,Vh) + stretch['U'] = np.dot(np.linalg.inv(R),F) + stretch['V'] = np.dot(F,np.linalg.inv(R)) + for theStretch in stretches: + for i in range(9): + if abs(stretch[theStretch][i%3,i//3]) < 1e-12: # kill nasty noisy data + stretch[theStretch][i%3,i//3] = 0.0 + (D,V) = np.linalg.eig(stretch[theStretch]) # eigen decomposition (of symmetric matrix) + for i,eigval in enumerate(D): + if eigval < 0.0: # flip negative eigenvalues + D[i] = -D[i] + V[:,i] = -V[:,i] + if np.dot(V[:,i],V[:,(i+1)%3]) != 0.0: # check each vector for orthogonality + V[:,(i+1)%3] = np.cross(V[:,(i+2)%3],V[:,i]) # correct next vector + V[:,(i+1)%3] /= np.sqrt(np.dot(V[:,(i+1)%3],V[:,(i+1)%3].conj())) # and renormalize (hyperphobic?) + for theStrain in strains: + d = operator(theStretch,theStrain,D) # operate on eigenvalues of U or V + eps = (np.dot(V,np.dot(np.diag(d),V.T)).real).reshape(9) # build tensor back from eigenvalue/vector basis - 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 norms - F = numpy.array(map(float,table.data[column['defgrad'][active['defgrad'][0]]: - column['defgrad'][active['defgrad'][0]]+datainfo['defgrad']['len']]),'d').reshape(3,3) - (U,S,Vh) = numpy.linalg.svd(F) - R = numpy.dot(U,Vh) - stretch['U'] = numpy.dot(numpy.linalg.inv(R),F) - stretch['V'] = numpy.dot(F,numpy.linalg.inv(R)) - for theStretch in stretches: - for i in range(9): - if abs(stretch[theStretch][i%3,i//3]) < 1e-12: # kill nasty noisy data - stretch[theStretch][i%3,i//3] = 0.0 - (D,V) = numpy.linalg.eig(stretch[theStretch]) # eigen decomposition (of symmetric matrix) - for i,eigval in enumerate(D): - if eigval < 0.0: # flip negative eigenvalues - D[i] = -D[i] - V[:,i] = -V[:,i] - if numpy.dot(V[:,i],V[:,(i+1)%3]) != 0.0: # check each vector for orthogonality - V[:,(i+1)%3] = numpy.cross(V[:,(i+2)%3],V[:,i]) # correct next vector - V[:,(i+1)%3] /= numpy.sqrt(numpy.dot(V[:,(i+1)%3],V[:,(i+1)%3].conj())) # and renormalize (hyperphobic?) - for theStrain in strains: - d = operator(theStretch,theStrain,D) # operate on eigenvalues of U or V - eps = (numpy.dot(V,numpy.dot(numpy.diag(d),V.T)).real).reshape(9) # build tensor back from eigenvalue/vector basis - - table.data_append(list(eps)) - - table.data_write() # output processed line + table.data_append(list(eps)) + outputAlive = table.data_write() # output processed line # ------------------------------------------ output result --------------------------------------- + outputAlive and table.output_flush() # just in case of buffered ASCII table - table.output_flush() # just in case of buffered ASCII table - - file['input'].close() # close input ASCII table + file['input'].close() # close input ASCII table (works for stdin) + file['output'].close() # close output ASCII table (works for stdout) if file['name'] != 'STDIN': - file['output'].close # close output ASCII table - os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new + os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new