#!/usr/bin/env python # -*- coding: UTF-8 no BOM -*- import os,sys,string import numpy as np from collections import defaultdict from optparse import OptionParser import damask scriptID = string.replace('$Id$','\n','\\n') scriptName = os.path.splitext(scriptID.split()[1])[0] # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ Add column(s) containing eigenvalues and eigenvectors of requested tensor column(s). """, version = scriptID) parser.add_option('-t','--tensor', dest = 'tensor', action = 'extend', metavar = '', help = 'heading of columns containing tensor field values') (options,filenames) = parser.parse_args() if options.tensor == None: parser.error('no data column specified.') # --- loop over input files ------------------------------------------------------------------------- if filenames == []: filenames = ['STDIN'] for name in filenames: if not (name == 'STDIN' or os.path.exists(name)): continue table = damask.ASCIItable(name = name, outname = name+'_tmp', buffered = False) table.croak('\033[1m'+scriptName+'\033[0m'+(': '+name if name != 'STDIN' else '')) # ------------------------------------------ read header ------------------------------------------ table.head_read() # ------------------------------------------ sanity checks ---------------------------------------- items = { 'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'column': []}, } errors = [] remarks = [] for type, data in items.iteritems(): for what in data['labels']: dim = table.label_dimension(what) if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type)) else: items[type]['column'].append(table.label_index(what)) table.labels_append(['{}_eigval({})'.format(i+1,what) for i in xrange(3)]) # extend ASCII header with new labels table.labels_append(['{}_eigvec({})'.format(i+1,what) for i in xrange(9)]) # extend ASCII header with new labels if remarks != []: table.croak(remarks) if errors != []: table.croak(errors) table.close(dismiss = True) continue # ------------------------------------------ assemble header -------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.head_write() # ------------------------------------------ process data ------------------------------------------ outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table for type, data in items.iteritems(): for column in data['column']: (u,v) = np.linalg.eig(np.array(map(float,table.data[column:column+data['dim']])).reshape(data['shape'])) table.data_append(list(u)) table.data_append(list(v.transpose().reshape(data['dim']))) outputAlive = table.data_write() # output processed line # ------------------------------------------ output finalization ----------------------------------- table.close() # close input ASCII table (works for stdin) if name != 'STDIN': os.rename(name+'_tmp',name) # overwrite old one with tmp new