From 43c18801954693e0130f27b708e48812725131f5 Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Tue, 29 Nov 2016 14:37:43 -0500
Subject: [PATCH] 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!
---
 processing/post/groupTable.py | 42 +++++++++++++++++------------------
 1 file changed, 20 insertions(+), 22 deletions(-)

diff --git a/processing/post/groupTable.py b/processing/post/groupTable.py
index a1205eff4..67d07a7d1 100755
--- a/processing/post/groupTable.py
+++ b/processing/post/groupTable.py
@@ -7,14 +7,11 @@ import numpy as np
 from optparse import OptionParser, OptionGroup
 import damask
 
-#"https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions"
-def periodicAverage(Points, Box):
-  theta = (Points/Box[1]) * (2.0*np.pi)
-  xi    = np.cos(theta)
-  zeta  = np.sin(theta)
-  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
+def periodicAverage(coords, limits):
+  """Centroid in periodic domain, see https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions"""
+  theta     = 2.0*np.pi * (coords - limits[0])/(limits[1] - limits[0])
+  theta_avg = np.pi + np.arctan2(-np.sin(theta).mean(axis=0), -np.cos(theta).mean(axis=0))
+  return limits[0] + theta_avg * (limits[1] - limits[0])/2.0/np.pi
 
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 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 = """
 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.
+Periodic domain averaging of coordinate values is supported.
 
 Examples:
 For grain averaged values, replace all rows of particular 'texture' with a single row containing their average.
@@ -43,23 +41,22 @@ parser.add_option('-a','--all',
                   dest = 'all',
                   action = 'store_true',
                   help = 'apply mapping function also to grouping column')
-
+                  
 group = OptionGroup(parser, "periodic averaging", "")
 
-group.add_option('-p','--periodic',
-                  dest = 'periodic',
-                  action = 'store_true',
-                  help = 'calculate average in periodic space defined by periodic length [%default]')
-group.add_option('--boundary',
-                 dest = 'boundary', metavar = 'MIN MAX',
-                 type = 'float', nargs = 2,
-                 help = 'define periodic box end points %default')
+group.add_option ('-p','--periodic',
+                  dest   = 'periodic',
+                  action = 'extend', metavar = '<string LIST>',
+                  help   = 'coordinate label(s) to average across periodic domain')
+group.add_option ('--limits',
+                  dest = 'boundary',
+                  type = 'float', metavar = 'float float', nargs = 2,
+                  help = 'min and max of periodic domain %default')
 
 parser.add_option_group(group)
 
 parser.set_defaults(function = 'np.average',
                     all      = False,
-                    periodic = False,
                     boundary = [0.0, 1.0])
 
 (options,filenames) = parser.parse_args()
@@ -108,6 +105,7 @@ for name in filenames:
 # ------------------------------------------ process data -------------------------------- 
 
   table.data_readArray()
+  indexrange = table.label_indexrange(options.periodic) if options.periodic is not None else []
   rows,cols  = table.data.shape
 
   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
   
   for i in range(len(values)):                                                                      # iterate over groups (unique values in grpColumn)
-    if options.periodic :
-       grpTable[i] = periodicAverage(table.data[index[i]:index[i+1]],options.boundary)              # apply periodicAverage mapping function
-    else :
-       grpTable[i] = np.apply_along_axis(mapFunction,0,table.data[index[i]:index[i+1]])             # apply mapping function
+    grpTable[i] = np.apply_along_axis(mapFunction,0,table.data[index[i]:index[i+1]])                # apply (general) mapping function
+    grpTable[i,indexrange] = \
+      periodicAverage(table.data[index[i]:index[i+1],indexrange],options.boundary)                  # apply periodicAverage mapping function
+
     if not options.all: grpTable[i,grpColumn] = table.data[index[i],grpColumn]                      # restore grouping column value
   
   table.data = grpTable