updated option --label for addCurl/Div/Grad

PRIVATE

@@ -1 +1 @@
-Subproject commit 1c1e8008489d81773a13a247604144a8d7ee3723
+Subproject commit bdbf2da71cd9e0825d17f673ec2fbabc2c8027f8

lib/damask/test/test.py

@@ -109,7 +109,7 @@ class Test():
 
         except Exception as e:
           logging.critical('exception during variant execution: "{}"'.format(str(e)))
-          return variant+1                                                     # return culprit
+          return variant+1                                                                            # return culprit
     return 0
   
   def feasible(self):

processing/post/addCurl.py

@@ -66,7 +66,7 @@ parser.add_option('-p','--pos','--periodiccellcenter',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-d','--data',
+parser.add_option('-l','--label',
                   dest = 'data',
                   action = 'extend', metavar = '<string LIST>',
                   help = 'label(s) of field values')

processing/post/addDivergence.py

@@ -9,36 +9,43 @@ import damask
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
+def merge_dicts(*dict_args):
+  """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts."""
+  result = {}
+  for dictionary in dict_args:
+      result.update(dictionary)
+  return result
+
 def divFFT(geomdim,field):
- shapeFFT    = np.array(np.shape(field))[0:3]
+  """Calculate divergence of a vector or tensor field by transforming into Fourier space."""
- grid = np.array(np.shape(field)[2::-1])
+  shapeFFT    = np.array(np.shape(field))[0:3]
- N = grid.prod()                                                                                    # field size
+  grid = np.array(np.shape(field)[2::-1])
- n = np.array(np.shape(field)[3:]).prod()                                                           # data size
+  N = grid.prod()                                                                                    # field size
+  n = np.array(np.shape(field)[3:]).prod()                                                           # data size
 
- if   n == 3:   dataType = 'vector'
+  field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
- elif n == 9:   dataType = 'tensor'
+  div_fourier   = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
 
- field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
+  # differentiation in Fourier space
- div_fourier   = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16')
+  TWOPIIMG = 2.0j*math.pi
+  einsums = { 
+              3:'ijkl,ijkl->ijk',                                                               # vector, 3 -> 1
+              9:'ijkm,ijklm->ijkl',                                                             # tensor, 3x3 -> 3
+            }
+  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                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
 
-# differentiation in Fourier space
+  k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
- TWOPIIMG = 2.0j*math.pi
+  if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # Nyquist freq=0 for even grid (Johnson, MIT, 2011)
- 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)
- 
- k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
- if grid[1]%2 == 0: k_sj[grid[1]//2] = 0                                                            # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
 
- k_si = np.arange(grid[0]//2+1)/geomdim[2]
+  k_si = np.arange(grid[0]//2+1)/geomdim[2]
- 
+
- kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
+  kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
- k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')                           
+  k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')                           
- if dataType == 'tensor':                                                                           # tensor, 3x3 -> 3
+
-   div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG
+  div_fourier = np.einsum(einsums[n],k_s,field_fourier)*TWOPIIMG
- elif dataType == 'vector':                                                                         # vector, 3 -> 1
+
-   div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG
+  return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
- 
- return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
 
 
 # --------------------------------------------------------------------
@@ -46,32 +53,38 @@ def divFFT(geomdim,field):
 # --------------------------------------------------------------------
 
 parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
-Add column(s) containing divergence of requested column(s).
+Add column(s) containing curl of requested column(s).
-Operates on periodic ordered three-dimensional data sets.
+Operates on periodic ordered three-dimensional data sets
-Deals with both vector- and tensor-valued fields.
+of vector and tensor fields.
-
 """, version = scriptID)
 
 parser.add_option('-p','--pos','--periodiccellcenter',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-v','--vector',
+parser.add_option('-l','--label',
-                  dest = 'vector',
+                  dest = 'data',
                   action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of vector field values')
+                  help = 'label(s) of field values')
-parser.add_option('-t','--tensor',
-                  dest = 'tensor',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'label(s) of tensor field values')
 
 parser.set_defaults(pos = 'pos',
                    )
 
+
 (options,filenames) = parser.parse_args()
 
-if options.vector is None and options.tensor is None:
+if options.data is None: parser.error('no data column specified.')
-  parser.error('no data column specified.')
+
+# --- define possible data types -------------------------------------------------------------------
+
+datatypes = {
+              3: {'name': 'vector',
+                  'shape': [3],
+                 },
+              9: {'name': 'tensor',
+                  'shape': [3,3],
+                 },
+            }
 
 # --- loop over input files ------------------------------------------------------------------------
 
@@ -82,30 +95,27 @@ for name in filenames:
   except: continue
   damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+# --- interpret header ----------------------------------------------------------------------------
 
   table.head_read()
 
-# ------------------------------------------ sanity checks ----------------------------------------
-
-  items = {
-            'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []},
-            'vector': {'dim': 3, 'shape': [3],   'labels':options.vector, 'active':[], 'column': []},
-          }
-  errors  = []
   remarks = []
-  column = {}
+  errors  = []
-  
+  active = []
-  if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
-  else: colCoord = table.label_index(options.pos)
 
-  for type, data in items.iteritems():
+  coordDim = table.label_dimension(options.pos)
-    for what in (data['labels'] if data['labels'] is not None else []):
+  if coordDim != 3:
-      dim = table.label_dimension(what)
+    errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos))
-      if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
+  else: coordCol = table.label_index(options.pos)
-      else:
+
-        items[type]['active'].append(what)
+  for me in options.data:
-        items[type]['column'].append(table.label_index(what))
+    dim = table.label_dimension(me)
+    if dim in datatypes:
+      active.append(merge_dicts({'label':me},datatypes[dim]))
+      remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me))
+    else:
+      remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \
+                     '"{}" not found...'.format(me) )
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:
@@ -116,17 +126,17 @@ for name in filenames:
 # ------------------------------------------ assemble header --------------------------------------
 
   table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  for type, data in items.iteritems():
+  for data in active:
-    for label in data['active']:
+    table.labels_append(['divFFT({})'.format(data['label']) if data['shape'] == [3] \
-      table.labels_append(['divFFT({})'.format(label) if type == 'vector' else
+                 else '{}_divFFT({})'.format(i+1,data['label']) 
-                           '{}_divFFT({})'.format(i+1,label) for i in range(data['dim']//3)])       # extend ASCII header with new labels
+                 for i in range(np.prod(np.array(data['shape']))//3)])                              # 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')
@@ -136,12 +146,11 @@ for name in filenames:
 # ------------------------------------------ process value field -----------------------------------
 
   stack = [table.data]
-  for type, data in items.iteritems():
+  for data in active:
-    for i,label in enumerate(data['active']):
+    # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
-      # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation
+    stack.append(divFFT(size[::-1],
-      stack.append(divFFT(size[::-1],
+                        table.data[:,table.label_indexrange(data['label'])].
-                          table.data[:,data['column'][i]:data['column'][i]+data['dim']].
+                        reshape(grid[::-1].tolist()+data['shape'])))
-                          reshape(grid[::-1].tolist()+data['shape'])))
 
 # ------------------------------------------ output result -----------------------------------------
 

processing/post/addGradient.py

@@ -63,7 +63,7 @@ parser.add_option('-p','--pos','--periodiccellcenter',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
                   help = 'label of coordinates [%default]')
-parser.add_option('-d','--data',
+parser.add_option('-l','--label',
                   dest = 'data',
                   action = 'extend', metavar = '<string LIST>',
                   help = 'label(s) of field values')