fixed periodic averaging to work with multi-dimensional data

option --periodic now takes list of labels that undergo periodoc domain averaging, i.e. incompatible to former API!
This commit is contained in:
Philip Eisenlohr 2016-11-29 14:37:43 -05:00
parent 7ebd688885
commit 43c1880195
1 changed files with 20 additions and 22 deletions

View File

@ -7,14 +7,11 @@ import numpy as np
from optparse import OptionParser, OptionGroup from optparse import OptionParser, OptionGroup
import damask import damask
#"https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions" def periodicAverage(coords, limits):
def periodicAverage(Points, Box): """Centroid in periodic domain, see https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions"""
theta = (Points/Box[1]) * (2.0*np.pi) theta = 2.0*np.pi * (coords - limits[0])/(limits[1] - limits[0])
xi = np.cos(theta) theta_avg = np.pi + np.arctan2(-np.sin(theta).mean(axis=0), -np.cos(theta).mean(axis=0))
zeta = np.sin(theta) return limits[0] + theta_avg * (limits[1] - limits[0])/2.0/np.pi
theta_avg = np.arctan2(-1.0*zeta.mean(), -1.0*xi.mean()) + np.pi
Pmean = Box[1] * theta_avg/(2.0*np.pi)
return Pmean
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
@ -26,6 +23,7 @@ scriptID = ' '.join([scriptName,damask.version])
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Apply a user-specified function to condense all rows for which column 'label' has identical values into a single row. Apply a user-specified function to condense all rows for which column 'label' has identical values into a single row.
Output table will contain as many rows as there are different (unique) values in the grouping column. Output table will contain as many rows as there are different (unique) values in the grouping column.
Periodic domain averaging of coordinate values is supported.
Examples: Examples:
For grain averaged values, replace all rows of particular 'texture' with a single row containing their average. For grain averaged values, replace all rows of particular 'texture' with a single row containing their average.
@ -46,20 +44,19 @@ parser.add_option('-a','--all',
group = OptionGroup(parser, "periodic averaging", "") group = OptionGroup(parser, "periodic averaging", "")
group.add_option('-p','--periodic', group.add_option ('-p','--periodic',
dest = 'periodic', dest = 'periodic',
action = 'store_true', action = 'extend', metavar = '<string LIST>',
help = 'calculate average in periodic space defined by periodic length [%default]') help = 'coordinate label(s) to average across periodic domain')
group.add_option('--boundary', group.add_option ('--limits',
dest = 'boundary', metavar = 'MIN MAX', dest = 'boundary',
type = 'float', nargs = 2, type = 'float', metavar = 'float float', nargs = 2,
help = 'define periodic box end points %default') help = 'min and max of periodic domain %default')
parser.add_option_group(group) parser.add_option_group(group)
parser.set_defaults(function = 'np.average', parser.set_defaults(function = 'np.average',
all = False, all = False,
periodic = False,
boundary = [0.0, 1.0]) boundary = [0.0, 1.0])
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
@ -108,6 +105,7 @@ for name in filenames:
# ------------------------------------------ process data -------------------------------- # ------------------------------------------ process data --------------------------------
table.data_readArray() table.data_readArray()
indexrange = table.label_indexrange(options.periodic) if options.periodic is not None else []
rows,cols = table.data.shape rows,cols = table.data.shape
table.data = table.data[np.lexsort([table.data[:,grpColumn]])] # sort data by grpColumn table.data = table.data[np.lexsort([table.data[:,grpColumn]])] # sort data by grpColumn
@ -117,10 +115,10 @@ for name in filenames:
grpTable = np.empty((len(values), cols)) # initialize output grpTable = np.empty((len(values), cols)) # initialize output
for i in range(len(values)): # iterate over groups (unique values in grpColumn) for i in range(len(values)): # iterate over groups (unique values in grpColumn)
if options.periodic : grpTable[i] = np.apply_along_axis(mapFunction,0,table.data[index[i]:index[i+1]]) # apply (general) mapping function
grpTable[i] = periodicAverage(table.data[index[i]:index[i+1]],options.boundary) # apply periodicAverage mapping function grpTable[i,indexrange] = \
else : periodicAverage(table.data[index[i]:index[i+1],indexrange],options.boundary) # apply periodicAverage mapping function
grpTable[i] = np.apply_along_axis(mapFunction,0,table.data[index[i]:index[i+1]]) # apply mapping function
if not options.all: grpTable[i,grpColumn] = table.data[index[i],grpColumn] # restore grouping column value if not options.all: grpTable[i,grpColumn] = table.data[index[i],grpColumn] # restore grouping column value
table.data = grpTable table.data = grpTable