diff --git a/PRIVATE b/PRIVATE index 2a19f3519..3341e5973 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2a19f35198e5e1e2f3e4d5a0728ed2667c2075f8 +Subproject commit 3341e5973bda63fe03ace6490dc6b010e188c3f3 diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index fefe8f84e..ef9d16504 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -4,6 +4,7 @@ import os,sys,math import numpy as np from optparse import OptionParser +from collections import defaultdict import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] @@ -22,6 +23,7 @@ def gradFFT(geomdim,field): grad_fourier = np.empty(field_fourier.shape+(3,),'c16') # differentiation in Fourier space +# Question: why are grid[0,1,2] normalized by geomdim[2,1,0]?? TWOPIIMG = 2.0j*math.pi k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0] if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) @@ -47,8 +49,8 @@ def gradFFT(geomdim,field): parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ Add column(s) containing gradient of requested column(s). -Operates on periodic ordered three-dimensional data sets. -Deals with both vector- and scalar fields. +Operates on periodic ordered three-dimensional data sets +of vector and scalar fields. """, version = scriptID) @@ -56,22 +58,17 @@ parser.add_option('-p','--pos','--periodiccellcenter', dest = 'pos', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') -parser.add_option('-v','--vector', - dest = 'vector', +parser.add_option('-d','--data', + dest = 'data', action = 'extend', metavar = '', - help = 'label(s) of vector field values') -parser.add_option('-s','--scalar', - dest = 'scalar', - action = 'extend', metavar = '', - help = 'label(s) of scalar field values') + help = 'label(s) of field values') parser.set_defaults(pos = 'pos', ) (options,filenames) = parser.parse_args() -if options.vector is None and options.scalar is None: - parser.error('no data column specified.') +if options.data is None: parser.error('no data column specified.') # --- loop over input files ------------------------------------------------------------------------ @@ -82,30 +79,31 @@ for name in filenames: except: continue damask.util.report(scriptName,name) -# ------------------------------------------ read header ------------------------------------------ +# --- interpret header ---------------------------------------------------------------------------- table.head_read() -# ------------------------------------------ sanity checks ---------------------------------------- - - items = { - 'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []}, - 'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []}, - } - errors = [] remarks = [] - column = {} - - if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos)) - else: colCoord = table.label_index(options.pos) + errors = [] + active = defaultdict(list) - for type, data in items.iteritems(): - for what in (data['labels'] if data['labels'] is not None else []): - dim = table.label_dimension(what) - if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) - else: - items[type]['active'].append(what) - items[type]['column'].append(table.label_index(what)) + coordDim = table.label_dimension(options.pos) + if coordDim != 3: + errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos)) + else: coordCol = table.label_index(options.pos) + + for i,dim in enumerate(table.label_dimension(options.data)): + me = options.data[i] + if dim == -1: + remarks.append('"{}" not found...'.format(me)) + elif dim == 1: + active['scalar'].append(me) + remarks.append('adding scalar "{}"...'.format(me)) + elif dim == 3: + active['vector'].append(me) + remarks.append('adding vector "{}"...'.format(me)) + else: + remarks.append('skipping "{}" of dimension {}...'.format(me,dim)) if remarks != []: damask.util.croak(remarks) if errors != []: @@ -113,19 +111,20 @@ for name in filenames: table.close(dismiss = True) continue + # ------------------------------------------ assemble header -------------------------------------- table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - for type, data in items.iteritems(): - for label in data['active']: - table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in range(3 * data['dim'])]) # extend ASCII header with new labels + for type, data in active.iteritems(): + for label in data: + table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in range(3*table.label_dimension(label))]) # extend ASCII header with new labels table.head_write() # --------------- figure out size and grid --------------------------------------------------------- table.data_readArray() - coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)] + coords = [np.unique(table.data[:,coordCol+i]) for i in range(3)] mincorner = np.array(map(min,coords)) maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') @@ -135,12 +134,12 @@ for name in filenames: # ------------------------------------------ process value field ----------------------------------- stack = [table.data] - for type, data in items.iteritems(): - for i,label in enumerate(data['active']): + for type, data in active.iteritems(): + for i,label in enumerate(data): # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation stack.append(gradFFT(size[::-1], - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. - reshape(grid[::-1].tolist()+data['shape']))) + table.data[:,table.label_indexrange(label)]. + reshape(grid[::-1].tolist()+[table.label_dimension(label)]))) # ------------------------------------------ output result -----------------------------------------