using central functionality

processing/post/addDisplacement.py

@@ -2,10 +2,10 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import numpy as np
-import scipy.ndimage
 
 import damask
 
@@ -14,79 +14,6 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
 
-#--------------------------------------------------------------------------------------------------
-def cell2node(cellData,grid):
-
-  nodeData = 0.0
-  datalen = np.array(cellData.shape[3:]).prod()
-  
-  for i in range(datalen):
-    node = scipy.ndimage.convolve(cellData.reshape(tuple(grid[::-1])+(datalen,))[...,i],
-                                  np.ones((2,2,2))/8.,                                              # 2x2x2 neighborhood of cells
-                                  mode = 'wrap',
-                                  origin = -1,                                                      # offset to have cell origin as center
-                                 )                                                                  # now averaged at cell origins
-    node = np.append(node,node[np.newaxis,0,:,:,...],axis=0)                                        # wrap along z
-    node = np.append(node,node[:,0,np.newaxis,:,...],axis=1)                                        # wrap along y
-    node = np.append(node,node[:,:,0,np.newaxis,...],axis=2)                                        # wrap along x
-
-    nodeData = node[...,np.newaxis] if i==0 else np.concatenate((nodeData,node[...,np.newaxis]),axis=-1)
-
-  return nodeData
-
-#--------------------------------------------------------------------------------------------------
-def displacementAvgFFT(F,grid,size,nodal=False,transformed=False):
-  """Calculate average cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
-  if nodal:
-    x, y, z = np.meshgrid(np.linspace(0,size[2],1+grid[2]),
-                          np.linspace(0,size[1],1+grid[1]),
-                          np.linspace(0,size[0],1+grid[0]),
-                          indexing = 'ij')
-  else:
-    delta = size/grid*0.5
-    x, y, z = np.meshgrid(np.linspace(delta[2],size[2]-delta[2],grid[2]),
-                          np.linspace(delta[1],size[1]-delta[1],grid[1]),
-                          np.linspace(delta[0],size[0]-delta[0],grid[0]),
-                          indexing = 'ij')
-
-  origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) 
-
-  F_fourier = F if transformed else np.fft.rfftn(F,axes=(0,1,2))                                    # transform or use provided data
-  Favg = np.real(F_fourier[0,0,0,:,:])/grid.prod()                                                  # take zero freq for average
-  avgDisplacement = np.einsum('ml,ijkl->ijkm',Favg-np.eye(3),origCoords)                            # dX = Favg.X
-
-  return avgDisplacement
-
-#--------------------------------------------------------------------------------------------------
-def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
-  """Calculate cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
-  integrator = 0.5j * size / np.pi
-
-  kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])),
-                           np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])),
-                                    np.arange(grid[0]//2+1),
-                           indexing = 'ij')
-  k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) 
-  k_sSquared = np.einsum('...l,...l',k_s,k_s)
-  k_sSquared[0,0,0] = 1.0                                                                           # ignore global average frequency
-
-#--------------------------------------------------------------------------------------------------
-# integration in Fourier space
-
-  displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm',
-                                    F if transformed else np.fft.rfftn(F,axes=(0,1,2)),
-                                    k_s,
-                                    integrator,
-                                   ) / k_sSquared[...,np.newaxis]
-
-#--------------------------------------------------------------------------------------------------
-# backtransformation to real space
-
-  displacement = np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2))
-
-  return cell2node(displacement,grid) if nodal else displacement
-
-
 # --------------------------------------------------------------------
 #                                MAIN
 # --------------------------------------------------------------------
@@ -100,7 +27,7 @@ Outputs at cell centers or cell nodes (into separate file).
 
 parser.add_option('-f',
                   '--defgrad',
-                  dest    = 'defgrad',
+                  dest    = 'f',
                   metavar = 'string',
                   help    = 'label of deformation gradient [%default]')
 parser.add_option('-p',
@@ -113,108 +40,35 @@ parser.add_option('--nodal',
                   action  = 'store_true',
                   help    = 'output nodal (instead of cell-centered) displacements')
 
-parser.set_defaults(defgrad = 'f',
+parser.set_defaults(f   = 'f',
                     pos = 'pos',
                    )
 
 (options,filenames) = parser.parse_args()
 
-# --- loop over input files -------------------------------------------------------------------------
-
-if filenames == []: filenames = [None]
-
 for name in filenames:
-  outname = (os.path.splitext(name)[0] +
+    damask.util.report(scriptName,name)
-             '_nodal' +
-             os.path.splitext(name)[1]) if (options.nodal and name) else None
-  try:    table = damask.ASCIItable(name = name,
-                                    outname = outname,
-                                    buffered = False)
-  except: continue
-  damask.util.report(scriptName,'{}{}'.format(name if name else '',
-                                              ' --> {}'.format(outname) if outname else ''))
 
-# ------------------------------------------ read header ------------------------------------------
+    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
+    grid,size = damask.util.coordGridAndSize(table.get_array(options.pos))
-  table.head_read()
-
-# ------------------------------------------ sanity checks ----------------------------------------
-
-  errors  = []
-  remarks = []
-  
-  if table.label_dimension(options.defgrad) != 9:
-    errors.append('deformation gradient "{}" is not a 3x3 tensor.'.format(options.defgrad))
-
-  coordDim = table.label_dimension(options.pos)
-  if not 3 >= coordDim >= 1:
-    errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
-  elif coordDim < 3:
-    remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim,
-                                                                            's' if coordDim < 2 else '',
-                                                                            options.pos))
-
-  if remarks != []: damask.util.croak(remarks)
-  if errors  != []:
-    damask.util.croak(errors)
-    table.close(dismiss=True)
-    continue
-
-# --------------- figure out size and grid ---------------------------------------------------------
-
-  table.data_readArray([options.defgrad,options.pos])
-  table.data_rewind()
-
-  if len(table.data.shape) < 2: table.data.shape += (1,)                                            # expand to 2D shape
-  if table.data[:,9:].shape[1] < 3:
-    table.data = np.hstack((table.data,
-                            np.zeros((table.data.shape[0],
-                                      3-table.data[:,9:].shape[1]),dtype='f')))                     # fill coords up to 3D with zeros
-
-  grid,size = damask.util.coordGridAndSize(table.data[:,9:12])
-  N = grid.prod()
-
-  if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid))
-  if errors  != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)
-    continue
-  
-# ------------------------------------------ process data ------------------------------------------
-
-  F_fourier = np.fft.rfftn(table.data[:,:9].reshape(grid[2],grid[1],grid[0],3,3),axes=(0,1,2))      # perform transform only once...
-
-  fluctDisplacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True)
-  avgDisplacement   = displacementAvgFFT  (F_fourier,grid,size,options.nodal,transformed=True)
-
-# ------------------------------------------ assemble header ---------------------------------------
     
+    F = table.get_array(options.f).reshape(np.append(grid[::-1],(3,3)))
     if options.nodal:
-    table.info_clear()
+        table = damask.Table(damask.grid_filters.coord0_node(grid[::-1],size[::-1]).reshape((-1,3)),
-    table.labels_clear()
+                             {'pos':(3,)})
+        table.add_array('avg({}).{}'.format(options.f,options.pos),
+                        damask.grid_filters.displacement_avg_node(size[::-1],F).reshape((-1,3)),
+                        scriptID+' '+' '.join(sys.argv[1:]))
+        table.add_array('fluct({}).{}'.format(options.f,options.pos),
+                        damask.grid_filters.displacement_fluct_node(size[::-1],F).reshape((-1,3)),
+                        scriptID+' '+' '.join(sys.argv[1:]))
+        table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt')
+    else:
+        table.add_array('avg({}).{}'.format(options.f,options.pos),
+                        damask.grid_filters.displacement_avg_cell(size[::-1],F).reshape((-1,3)),
+                        scriptID+' '+' '.join(sys.argv[1:]))
+        table.add_array('fluct({}).{}'.format(options.f,options.pos),
+                        damask.grid_filters.displacement_fluct_cell(size[::-1],F).reshape((-1,3)),
+                        scriptID+' '+' '.join(sys.argv[1:]))
     
-  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
+        table.to_ASCII(sys.stdout if name is None else name)
-  table.labels_append((['{}_pos'         .format(i+1)   for i in range(3)] if options.nodal else []) +
-                       ['{}_avg({}).{}'  .format(i+1,options.defgrad,options.pos) for i in range(3)] +
-                       ['{}_fluct({}).{}'.format(i+1,options.defgrad,options.pos) for i in range(3)] )
-  table.head_write()
-
-# ------------------------------------------ output data -------------------------------------------
-
-  Zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else range(grid[2])
-  Yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else range(grid[1])
-  Xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else range(grid[0])
-
-  for i,z     in enumerate(Zrange):
-    for j,y   in enumerate(Yrange):
-      for k,x in enumerate(Xrange):
-        if options.nodal: table.data_clear()
-        else:             table.data_read()
-        table.data_append([x,y,z] if options.nodal else [])
-        table.data_append(list(  avgDisplacement[i,j,k,:]))
-        table.data_append(list(fluctDisplacement[i,j,k,:]))
-        table.data_write()                       
-
-# ------------------------------------------ output finalization -----------------------------------
-
-  table.close()                                                                                     # close ASCII tables

python/damask/grid_filters.py

@@ -102,14 +102,21 @@ def displacement_fluct_node(size,F):
 
 def displacement_avg_node(size,F):
     F_avg = np.average(F,axis=(0,1,2))
-    return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),coord0_node(F.shape[0:3],size))
+    return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),coord0_node(F.shape[:3],size))
 
 
 def cell_2_node(cell_data):
     """Interpolate cell data to nodal data."""
-  
     n = (  cell_data + np.roll(cell_data,1,(0,1,2))
          + np.roll(cell_data,1,(0,))  + np.roll(cell_data,1,(1,))  + np.roll(cell_data,1,(2,))
-         + np.roll(cell_data,1,(0,1)) + np.roll(cell_data,1,(1,2)) + np.roll(cell_data,1,(2,0))) *0.125
+         + np.roll(cell_data,1,(0,1)) + np.roll(cell_data,1,(1,2)) + np.roll(cell_data,1,(2,0)))*0.125
   
     return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
+
+def node_2_cell(node_data):
+    """Interpolate nodal data to cell data."""
+    c = (  node_data + np.roll(node_data,1,(0,1,2))
+         + np.roll(node_data,1,(0,))  + np.roll(node_data,1,(1,))  + np.roll(node_data,1,(2,))
+         + np.roll(node_data,1,(0,1)) + np.roll(node_data,1,(1,2)) + np.roll(node_data,1,(2,0)))*0.125
+  
+    return c[:-1,:-1,:-1]