From f49091c952c00f561337871fa7b21943766e08c4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 2 Dec 2015 20:53:02 +0000 Subject: [PATCH] addDeformedConfiguration without core module --- processing/post/addDeformedConfiguration.py | 103 ++++++++++++++------ 1 file changed, 71 insertions(+), 32 deletions(-) diff --git a/processing/post/addDeformedConfiguration.py b/processing/post/addDeformedConfiguration.py index 167ec2c56..a83408334 100755 --- a/processing/post/addDeformedConfiguration.py +++ b/processing/post/addDeformedConfiguration.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: UTF-8 no BOM -*- -import os,sys,string +import os,sys,string,math import numpy as np from optparse import OptionParser import damask @@ -9,6 +9,49 @@ import damask scriptID = string.replace('$Id$','\n','\\n') scriptName = os.path.splitext(scriptID.split()[1])[0] +#-------------------------------------------------------------------------------------------------- +def deformedCoordsFFT(F,undeformed=False): +#-------------------------------------------------------------------------------------------------- + wgt = 1.0/grid.prod() + integrator = grid[::-1] / 2.0 / math.pi + step = size[::-1]/grid[::-1] + + F_fourier = np.fft.fftpack.rfftn(F,axes=(0,1,2)) + if undeformed: + Favg=np.eye(3) + else: + Favg=np.real(F_fourier[0,0,0,:,:])*wgt + coords_fourier = np.zeros(F_fourier.shape[:-1],'c16') + +#-------------------------------------------------------------------------------------------------- +# integration in Fourier space + k_s = np.zeros([3],'i') + + for k in xrange(grid[2]): + k_s[2] = k + if(k > grid[0]/2 ): k_s[2] = k_s[2] - grid[2] + for j in xrange(grid[1]): + k_s[1] = j + if(j > grid[1]/2 ): k_s[1] = k_s[1] - grid[1] + for i in xrange(grid[0]/2+1): + k_s[0] = i + for m in xrange(3): + coords_fourier[k,j,i,m] = sum(F_fourier[k,j,i,0:3,m]*k_s*\ + np.array([0.+1.j,0.+1.j,0.+1.j],'c16')*integrator) + if (any(k_s != 0)): + coords_fourier[k,j,i,0:3] /= -sum(k_s*k_s) +#-------------------------------------------------------------------------------------------------- +# add average to scaled fluctuation and put (0,0,0) on (0,0,0) + coords = np.fft.fftpack.irfftn(coords_fourier,axes=(0,1,2)) + offset_coords = np.dot(F[0,0,0,0:3,0:3],step/2.0) - scaling*coords[0,0,0,0:3] + for k in xrange(grid[2]): + for j in xrange(grid[1]): + for i in xrange(grid[0]): + coords[k,j,i,0:3] = scaling*coords[k,j,i,0:3] \ + + offset_coords \ + + np.dot(Favg,step*np.array([i,j,k])) + return coords + # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -19,26 +62,29 @@ Operates on periodic three-dimensional x,y,z-ordered data sets. """, version = scriptID) -parser.add_option('-c','--coordinates', - dest = 'coords', - type = 'string', metavar = 'string', - help = 'column label of coordinates [%default]') -parser.add_option('-f','--defgrad', - dest = 'defgrad', - type = 'string', metavar = 'string', - help = 'column label of deformation gradient [%default]') -parser.add_option('--scaling', - dest = 'scaling', - type = 'float', nargs = 3, metavar = ' '.join(['float']*3), - help = 'x/y/z scaling of displacement fluctuation') +parser.add_option('-f', '--defgrad',dest='defgrad', metavar = 'string', + help='heading of deformation gradient columns [%default]') +parser.add_option('--reference', dest='undeformed', action='store_true', + help='map results to reference (undeformed) average configuration [%default]') +parser.add_option('--scaling', dest='scaling', action='extend', metavar = '', + help='scaling of fluctuation') +parser.add_option('-u', '--unitlength', dest='unitlength', type='float', metavar = 'float', + help='set unit length for 2D model [%default]') +parser.add_option('--coordinates', dest='coords', metavar='string', + help='column heading for coordinates [%default]') -parser.set_defaults(coords = 'ipinitialcoord', - defgrad = 'f', - scaling = [1.,1.,1.], - ) +parser.set_defaults(defgrad = 'f') +parser.set_defaults(coords = 'ipinitialcoord') +parser.set_defaults(scaling = []) +parser.set_defaults(undeformed = False) +parser.set_defaults(unitlength = 0.0) (options,filenames) = parser.parse_args() +options.scaling += [1.0 for i in xrange(max(0,3-len(options.scaling)))] +scaling = map(float, options.scaling) + + # --- loop over input files ------------------------------------------------------------------------- if filenames == []: filenames = [None] @@ -83,7 +129,7 @@ for name in filenames: size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other spacings 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) @@ -94,23 +140,16 @@ for name in filenames: table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) for coord in xrange(3): - table.labels_append(['{}_{}.{}'.format(coord+1,options.defgrad,options.coords) ]) # extend ASCII header with new labels + label = '{}_{}_{}'.format(coord+1,options.defgrad,options.coords) + if np.any(scaling) != 1.0: label+='_{}_{}_{}'.format(scaling) + if options.undeformed: label+='_undeformed' + table.labels_append([label]) # extend ASCII header with new labels table.head_write() # ------------------------------------------ read deformation gradient field ----------------------- - table.data_rewind() - F = np.array([0.0 for i in xrange(N*9)]).reshape([3,3]+grid.tolist()) - idx = 0 - while table.data_read(): - (x,y,z) = damask.util.gridLocation(idx,grid) # figure out (x,y,z) position from line count - idx += 1 - F[0:3,0:3,x,y,z] = np.array(map(float,table.data[table.label_index(options.defgrad):\ - table.label_index(options.defgrad)+9]),'d').reshape(3,3) + centroids = deformedCoordsFFT(table.data[:,colF:colF+9].reshape(grid[0],grid[1],grid[2],3,3), + options.undeformed) -# ------------------------------------------ calculate coordinates --------------------------------- - Favg = damask.core.math.tensorAvg(F) - centroids = damask.core.mesh.deformedCoordsFFT(size,F,Favg) - # ------------------------------------------ process data ------------------------------------------ table.data_rewind() idx = 0 @@ -118,7 +157,7 @@ for name in filenames: while outputAlive and table.data_read(): # read next data line of ASCII table (x,y,z) = damask.util.gridLocation(idx,grid) # figure out (x,y,z) position from line count idx += 1 - table.data_append(list(centroids[:,x,y,z])) + table.data_append(list(centroids[z,y,x,:])) outputAlive = table.data_write() # ------------------------------------------ output finalization -----------------------------------