From 69c1cb71de23f0fb1fc5b616fddedc05c3e11c05 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 9 Mar 2016 04:12:09 +0100 Subject: [PATCH 01/18] updated version information after successful test of v2.0.0-6-gb5302ec --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 2bc4b2cfc..d9e8f873a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-3-g618e352 +v2.0.0-6-gb5302ec From bb9ee3910fac62010ab646cf3aaa841d0ccf6d7e Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 9 Mar 2016 16:12:54 +0100 Subject: [PATCH 02/18] updated version information after successful test of v2.0.0-7-g2296929 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 2bc4b2cfc..4c6293bf0 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-3-g618e352 +v2.0.0-7-g2296929 From 2dbdf9a72f28f488746b598b9914ae008b03b455 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 13 Mar 2016 18:21:25 +0100 Subject: [PATCH 03/18] all symbols need to be representable in unicode --- lib/damask/util.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/damask/util.py b/lib/damask/util.py index 77bc30403..9ce42759a 100644 --- a/lib/damask/util.py +++ b/lib/damask/util.py @@ -104,7 +104,8 @@ class extendableOption(Option): class backgroundMessage(threading.Thread): """reporting with animation to indicate progress""" - choices = {'bounce': ['_','o','O','°','¯','¯','°','O','o',], + choices = {'bounce': ['_', 'o', 'O', u'\u00B0', + u'\u203e',u'\u203e',u'\u00B0','O','o','_'], 'circle': [u'\u25f4',u'\u25f5',u'\u25f6',u'\u25f7'], 'hexagon': [u'\u2b22',u'\u2b23'], 'square': [u'\u2596',u'\u2598',u'\u259d',u'\u2597'], From 80d2bdbfde1c55fbf3332eeedfd8565efa19de65 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 14 Mar 2016 04:12:59 +0100 Subject: [PATCH 04/18] updated version information after successful test of v2.0.0-18-g2dbdf9a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index c02195acf..6a337695e 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-16-g344c6f6 +v2.0.0-18-g2dbdf9a From b40e7a26124a59b8d38240d4b49a2012e0a42f3d Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 14 Mar 2016 16:11:55 +0100 Subject: [PATCH 05/18] updated version information after successful test of v2.0.0-19-g80d2bdb --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 6a337695e..9f18a71da 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-18-g2dbdf9a +v2.0.0-19-g80d2bdb From 6197e839c942b5dc91fb7243bc0423fff39b3264 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 15 Mar 2016 14:52:10 -0400 Subject: [PATCH 06/18] fixed missing file extension when adding in-place --- processing/post/vtk_addRectilinearGridData.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/vtk_addRectilinearGridData.py b/processing/post/vtk_addRectilinearGridData.py index 922b4b696..3ca54d84e 100755 --- a/processing/post/vtk_addRectilinearGridData.py +++ b/processing/post/vtk_addRectilinearGridData.py @@ -151,7 +151,7 @@ for name in filenames: writer = vtk.vtkXMLRectilinearGridWriter() writer.SetDataModeToBinary() writer.SetCompressorTypeToZLib() - writer.SetFileName(os.path.splitext(options.vtk)[0]+('' if options.inplace else '_added.vtr')) + writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtr' if options.inplace else '_added.vtr')) if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid) else: writer.SetInputData(rGrid) writer.Write() From 9704eb6926f6be905128ae0552dec146fa811eb7 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 15 Mar 2016 17:23:01 -0400 Subject: [PATCH 07/18] small correction to variable declaration layout --- code/lattice.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/code/lattice.f90 b/code/lattice.f90 index f5717e8e1..05a123125 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -2182,7 +2182,7 @@ pure function lattice_qDisorientation(Q1, Q2, struct) real(pReal), dimension(4) :: lattice_qDisorientation real(pReal), dimension(4), intent(in) :: & Q1, & ! 1st orientation - Q2 ! 2nd orientation + Q2 ! 2nd orientation integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered struct From 93466273ffd03db9cc03e6fadf02727ba3967d19 Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 16 Mar 2016 04:13:16 +0100 Subject: [PATCH 08/18] updated version information after successful test of v2.0.0-31-g6197e83 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 9f18a71da..27bd816bb 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-19-g80d2bdb +v2.0.0-31-g6197e83 From ce26ad06da8cbc945444412d42220b4755562730 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Wed, 16 Mar 2016 14:55:56 -0400 Subject: [PATCH 09/18] polishing and introduction of locally derived grid in FFT subroutine --- processing/post/addCurl.py | 39 +++++++++++++------------------- processing/post/addDivergence.py | 36 ++++++++++++----------------- 2 files changed, 31 insertions(+), 44 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index ac9715c36..6789b27d5 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -10,40 +10,35 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def curlFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) N = grid.prod() # field size n = np.array(np.shape(field)[3:]).prod() # data size - if n == 3: - dataType = 'vector' - elif n == 9: - dataType = 'tensor' + if n == 3: dataType = 'vector' + elif n == 9: dataType = 'tensor' field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) curl_fourier = np.zeros(field_fourier.shape,'c16') # differentiation in Fourier space k_s = np.zeros([3],'i') - TWOPIIMG = (0.0+2.0j*math.pi) + TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if(grid[2]%2==0 and i == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[0]=0 - elif (i > grid[2]//2): - k_s[0] = k_s[0] - grid[2] + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[1]=0 - elif (j > grid[1]//2): - k_s[1] = k_s[1] - grid[1] + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if(grid[0]%2==0 and k == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[2]=0 + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order - xi = np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16') if dataType == 'tensor': for l in xrange(3): curl_fourier[i,j,k,0,l] = ( field_fourier[i,j,k,l,2]*xi[1]\ @@ -100,10 +95,8 @@ if options.vector is None and options.tensor is None: if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name,buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ @@ -161,9 +154,9 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(curlFFT(size[::-1], # we need to reverse order here, because x - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. # is fastest,ie rightmost, but leftmost in - reshape([grid[2],grid[1],grid[0]]+data['shape']))) # our x,y,z notation + stack.append(curlFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) # ------------------------------------------ output result ----------------------------------------- diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 75708a603..99b56ead2 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -10,6 +10,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def divFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) N = grid.prod() # field size n = np.array(np.shape(field)[3:]).prod() # data size @@ -18,27 +19,22 @@ def divFFT(geomdim,field): # differentiation in Fourier space k_s=np.zeros([3],'i') - TWOPIIMG = (0.0+2.0j*math.pi) + TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if(grid[2]%2==0 and i == grid[2]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[0]=0 - elif (i > grid[2]//2): - k_s[0] = k_s[0] - grid[2] + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if(grid[1]%2==0 and j == grid[1]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[1]=0 - elif (j > grid[1]//2): - k_s[1] = k_s[1] - grid[1] + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if(grid[0]%2==0 and k == grid[0]//2): # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - k_s[2]=0 + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi=np.array([k_s[2]/geomdim[2]+0.0j,k_s[1]/geomdim[1]+0.j,k_s[0]/geomdim[0]+0.j],'c16') + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order if n == 9: # tensor, 3x3 -> 3 for l in xrange(3): div_fourier[i,j,k,l] = sum(field_fourier[i,j,k,l,0:3]*xi) *TWOPIIMG @@ -85,10 +81,8 @@ if options.vector is None and options.tensor is None: if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name,buffered = False) - except: - continue + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ @@ -140,16 +134,16 @@ for name in filenames: maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 set to smallest among other spacings + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other spacings # ------------------------------------------ process value field ----------------------------------- stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(divFFT(size[::-1], # we need to reverse order here, because x - table.data[:,data['column'][i]:data['column'][i]+data['dim']]. # is fastest,ie rightmost, but leftmost in - reshape([grid[2],grid[1],grid[0]]+data['shape']))) # our x,y,z notation + stack.append(divFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) # ------------------------------------------ output result ----------------------------------------- @@ -158,4 +152,4 @@ for name in filenames: # ------------------------------------------ output finalization ----------------------------------- - table.close() # close input ASCII table (works for stdin) \ No newline at end of file + table.close() # close input ASCII table (works for stdin) From 3a267dc7ffc4f65d18bd30f9743c171e9d389d41 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Wed, 16 Mar 2016 15:02:38 -0400 Subject: [PATCH 10/18] FFT-based gradient calculation of scalar or vector field data (sibling of addCurl and addDivergence) --- processing/post/addGradient.py | 157 +++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100755 processing/post/addGradient.py diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py new file mode 100755 index 000000000..0d7caa829 --- /dev/null +++ b/processing/post/addGradient.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python +# -*- coding: UTF-8 no BOM -*- + +import os,sys,math +import numpy as np +from optparse import OptionParser +import damask + +scriptName = os.path.splitext(os.path.basename(__file__))[0] +scriptID = ' '.join([scriptName,damask.version]) + +def gradFFT(geomdim,field): + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size + if n == 3: dataType = 'vector' + elif n == 1: dataType = 'scalar' + + field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) + grad_fourier = np.zeros(field_fourier.shape+(3,),'c16') + +# differentiation in Fourier space + k_s = np.zeros([3],'i') + TWOPIIMG = 2.0j*math.pi + for i in xrange(grid[2]): + k_s[0] = i + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif i > grid[2]//2: k_s[0] -= grid[2] + + for j in xrange(grid[1]): + k_s[1] = j + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + elif j > grid[1]//2: k_s[1] -= grid[1] + + for k in xrange(grid[0]//2+1): + k_s[2] = k + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order + + grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data + + if dataType == 'vector': + grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data + grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG + + return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2)).reshape([N,3*n]) + + +# -------------------------------------------------------------------- +# MAIN +# -------------------------------------------------------------------- + +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ +Add column(s) containing gradient of requested column(s). +Operates on periodic ordered three-dimensional data sets. +Deals with both vector- and scalar fields. + +""", version = scriptID) + +parser.add_option('-c','--coordinates', + dest = 'coords', + type = 'string', metavar='string', + help = 'column heading for coordinates [%default]') +parser.add_option('-v','--vector', + dest = 'vector', + action = 'extend', metavar = '', + help = 'heading of columns containing vector field values') +parser.add_option('-s','--scalar', + dest = 'scalar', + action = 'extend', metavar = '', + help = 'heading of columns containing scalar field values') + +parser.set_defaults(coords = 'ipinitialcoord', + ) + +(options,filenames) = parser.parse_args() + +if options.vector is None and options.scalar is None: + parser.error('no data column specified.') + +# --- loop over input files ------------------------------------------------------------------------- + +if filenames == []: filenames = [None] + +for name in filenames: + try: table = damask.ASCIItable(name = name,buffered = False) + except: continue + damask.util.report(scriptName,name) + +# ------------------------------------------ read header ------------------------------------------ + + table.head_read() + +# ------------------------------------------ sanity checks ---------------------------------------- + + items = { + 'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []}, + 'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []}, + } + errors = [] + remarks = [] + column = {} + + if table.label_dimension(options.coords) != 3: errors.append('coordinates {} are not a vector.'.format(options.coords)) + else: colCoord = table.label_index(options.coords) + + for type, data in items.iteritems(): + for what in (data['labels'] if data['labels'] is not None else []): + dim = table.label_dimension(what) + if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type)) + else: + items[type]['active'].append(what) + items[type]['column'].append(table.label_index(what)) + + if remarks != []: damask.util.croak(remarks) + if errors != []: + damask.util.croak(errors) + table.close(dismiss = True) + continue + +# ------------------------------------------ assemble header -------------------------------------- + + table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) + for type, data in items.iteritems(): + for label in data['active']: + table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels # grad increases the field dimension by one + table.head_write() + +# --------------- figure out size and grid --------------------------------------------------------- + + table.data_readArray() + + coords = [np.unique(table.data[:,colCoord+i]) for i in xrange(3)] + mincorner = np.array(map(min,coords)) + maxcorner = np.array(map(max,coords)) + grid = np.array(map(len,coords),'i') + size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) + +# ------------------------------------------ process value field ----------------------------------- + + stack = [table.data] + for type, data in items.iteritems(): + for i,label in enumerate(data['active']): + stack.append(gradFFT(size[::-1], # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation + table.data[:,data['column'][i]:data['column'][i]+data['dim']]. + reshape([grid[2],grid[1],grid[0]]+data['shape']))) + +# ------------------------------------------ output result ----------------------------------------- + + if len(stack) > 1: table.data = np.hstack(tuple(stack)) + table.data_writeArray('%.12g') + +# ------------------------------------------ output finalization ----------------------------------- + + table.close() # close input ASCII table (works for stdin) From 36c703ecd34f0c9033f372e6799dd831b6daede7 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Wed, 16 Mar 2016 15:12:53 -0400 Subject: [PATCH 11/18] shortened excessively long lines --- processing/post/addCurl.py | 15 ++++++++------- processing/post/addDivergence.py | 25 +++++++++++++------------ processing/post/addGradient.py | 23 ++++++++++++----------- 3 files changed, 33 insertions(+), 30 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 6789b27d5..d8b1ee025 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -11,8 +11,8 @@ scriptID = ' '.join([scriptName,damask.version]) def curlFFT(geomdim,field): grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - 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' elif n == 9: dataType = 'tensor' @@ -25,19 +25,19 @@ def curlFFT(geomdim,field): TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order if dataType == 'tensor': for l in xrange(3): @@ -154,7 +154,8 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(curlFFT(size[::-1], # 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(curlFFT(size[::-1], table.data[:,data['column'][i]:data['column'][i]+data['dim']]. reshape([grid[2],grid[1],grid[0]]+data['shape']))) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 99b56ead2..aadaceabf 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -11,34 +11,34 @@ scriptID = ' '.join([scriptName,damask.version]) def divFFT(geomdim,field): grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - 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 field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2)) - div_fourier = np.zeros(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector + div_fourier = np.zeros(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector # differentiation in Fourier space k_s=np.zeros([3],'i') TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order - if n == 9: # tensor, 3x3 -> 3 + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field input order + if n == 9: # tensor, 3x3 -> 3 for l in xrange(3): div_fourier[i,j,k,l] = sum(field_fourier[i,j,k,l,0:3]*xi) *TWOPIIMG - elif n == 3: # vector, 3 -> 1 + elif n == 3: # vector, 3 -> 1 div_fourier[i,j,k] = sum(field_fourier[i,j,k,0:3]*xi) *TWOPIIMG return np.fft.fftpack.irfftn(div_fourier,axes=(0,1,2)).reshape([N,n/3]) @@ -76,7 +76,7 @@ parser.set_defaults(coords = 'ipinitialcoord', if options.vector is None and options.tensor is None: parser.error('no data column specified.') -# --- loop over input files ------------------------------------------------------------------------- +# --- loop over input files ------------------------------------------------------------------------ if filenames == []: filenames = [None] @@ -134,14 +134,15 @@ for name in filenames: maxcorner = np.array(map(max,coords)) grid = np.array(map(len,coords),'i') size = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner) # size from edge to edge = dim * n/(n-1) - size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other spacings + size = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1])) # spacing for grid==1 equal to smallest among other ones # ------------------------------------------ process value field ----------------------------------- stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(divFFT(size[::-1], # 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], table.data[:,data['column'][i]:data['column'][i]+data['dim']]. reshape([grid[2],grid[1],grid[0]]+data['shape']))) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 0d7caa829..555e587de 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -11,8 +11,8 @@ scriptID = ' '.join([scriptName,damask.version]) def gradFFT(geomdim,field): grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - 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' elif n == 1: dataType = 'scalar' @@ -24,24 +24,24 @@ def gradFFT(geomdim,field): TWOPIIMG = 2.0j*math.pi for i in xrange(grid[2]): k_s[0] = i - if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[2]%2 == 0 and i == grid[2]//2: k_s[0] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif i > grid[2]//2: k_s[0] -= grid[2] for j in xrange(grid[1]): k_s[1] = j - if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[1]%2 == 0 and j == grid[1]//2: k_s[1] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) elif j > grid[1]//2: k_s[1] -= grid[1] for k in xrange(grid[0]//2+1): k_s[2] = k - if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) + if grid[0]%2 == 0 and k == grid[0]//2: k_s[2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011) - xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order + xi = (k_s/geomdim)[2::-1].astype('c16') # reversing the field order - grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data + grad_fourier[i,j,k,0,:] = field_fourier[i,j,k,0]*xi *TWOPIIMG # vector field from scalar data if dataType == 'vector': - grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data + grad_fourier[i,j,k,1,:] = field_fourier[i,j,k,1]*xi *TWOPIIMG # tensor field from vector data grad_fourier[i,j,k,2,:] = field_fourier[i,j,k,2]*xi *TWOPIIMG return np.fft.fftpack.irfftn(grad_fourier,axes=(0,1,2)).reshape([N,3*n]) @@ -79,7 +79,7 @@ parser.set_defaults(coords = 'ipinitialcoord', if options.vector is None and options.scalar is None: parser.error('no data column specified.') -# --- loop over input files ------------------------------------------------------------------------- +# --- loop over input files ------------------------------------------------------------------------ if filenames == []: filenames = [None] @@ -124,7 +124,7 @@ for name in filenames: table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) for type, data in items.iteritems(): for label in data['active']: - table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels # grad increases the field dimension by one + table.labels_append(['{}_gradFFT({})'.format(i+1,label) for i in xrange(3 * data['dim'])]) # extend ASCII header with new labels table.head_write() # --------------- figure out size and grid --------------------------------------------------------- @@ -143,7 +143,8 @@ for name in filenames: stack = [table.data] for type, data in items.iteritems(): for i,label in enumerate(data['active']): - stack.append(gradFFT(size[::-1], # 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(gradFFT(size[::-1], table.data[:,data['column'][i]:data['column'][i]+data['dim']]. reshape([grid[2],grid[1],grid[0]]+data['shape']))) From 36c88cde02a62a0e13961dcf0227907ee93e0a2f Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 17 Mar 2016 04:13:35 +0100 Subject: [PATCH 12/18] updated version information after successful test of v2.0.0-36-g2b7524e --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 27bd816bb..09b18f22c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-31-g6197e83 +v2.0.0-36-g2b7524e From d42f8d5fd4958bbc1015f4db495995e814fbf823 Mon Sep 17 00:00:00 2001 From: Zhuowen Zhao Date: Thu, 17 Mar 2016 17:40:40 -0400 Subject: [PATCH 13/18] added possibility to read real numbers as microstructure indices (not only integers) --- lib/damask/asciitable.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 4abe4297d..5212ef00c 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -517,21 +517,25 @@ class ASCIItable(): # ------------------------------------------------------------------ def microstructure_read(self, - grid): + grid, + type = 'i'): """read microstructure data (from .geom format)""" - N = grid.prod() # expected number of microstructure indices in data - microstructure = np.zeros(N,'i') # initialize as flat array + def datatype(item): + return int(item) if type.lower() == 'i' else float(item) + + N = grid.prod() # expected number of microstructure indices in data + microstructure = np.zeros(N,type) # initialize as flat array i = 0 while i < N and self.data_read(): items = self.data if len(items) > 2: - if items[1].lower() == 'of': items = [int(items[2])]*int(items[0]) - elif items[1].lower() == 'to': items = range(int(items[0]),1+int(items[2])) - else: items = map(int,items) - else: items = map(int,items) + if items[1].lower() == 'of': items = np.ones(datatype(items[0]))*datatype(items[2]) + elif items[1].lower() == 'to': items = np.arange(datatype(items[0]),1+datatype(items[2])) + else: items = map(datatype,items) + else: items = map(datatype,items) - s = min(len(items), N-i) # prevent overflow of microstructure array + s = min(len(items), N-i) # prevent overflow of microstructure array microstructure[i:i+s] = items[:s] i += s From 0f78f1c74713b389462e2aa745a158c2cdcbfd45 Mon Sep 17 00:00:00 2001 From: Aritra Chakraborty Date: Thu, 17 Mar 2016 18:09:22 -0400 Subject: [PATCH 14/18] possibility to treat floating point numbers as microstructure indices --- processing/pre/geom_canvas.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/processing/pre/geom_canvas.py b/processing/pre/geom_canvas.py index 882f32c88..ec171f522 100755 --- a/processing/pre/geom_canvas.py +++ b/processing/pre/geom_canvas.py @@ -28,24 +28,32 @@ parser.add_option('-o', '--offset', help = 'a,b,c offset from old to new origin of grid [%default]') parser.add_option('-f', '--fill', dest = 'fill', - type = 'int', metavar = 'int', + type = 'float', metavar = 'float', help = '(background) canvas grain index. "0" selects maximum microstructure index + 1 [%default]') +parser.add_option('--float', + dest = 'real', + action = 'store_true', + help = 'input data is float [%default]') parser.set_defaults(grid = ['0','0','0'], offset = (0,0,0), fill = 0, + real = False, ) (options, filenames) = parser.parse_args() +datatype = 'f' if options.real else 'i' + + # --- loop over input files ------------------------------------------------------------------------- if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name, - buffered = False, labeled = False) + try: table = damask.ASCIItable(name = name, + buffered = False, + labeled = False) except: continue damask.util.report(scriptName,name) @@ -71,7 +79,7 @@ for name in filenames: # --- read data ------------------------------------------------------------------------------------ - microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F') # read microstructure + microstructure = table.microstructure_read(info['grid'],datatype).reshape(info['grid'],order='F') # read microstructure # --- do work ------------------------------------------------------------------------------------ @@ -85,8 +93,8 @@ for name in filenames: else int(n) for o,n in zip(info['grid'],options.grid)],'i') newInfo['grid'] = np.where(newInfo['grid'] > 0, newInfo['grid'],info['grid']) - microstructure_cropped = np.zeros(newInfo['grid'],'i') - microstructure_cropped.fill(options.fill if options.fill > 0 else microstructure.max()+1) + microstructure_cropped = np.zeros(newInfo['grid'],datatype) + microstructure_cropped.fill(options.fill if options.real or options.fill > 0 else microstructure.max()+1) xindex = list(set(xrange(options.offset[0],options.offset[0]+newInfo['grid'][0])) & \ set(xrange(info['grid'][0]))) yindex = list(set(xrange(options.offset[1],options.offset[1]+newInfo['grid'][1])) & \ @@ -152,9 +160,9 @@ for name in filenames: # --- write microstructure information ------------------------------------------------------------ - formatwidth = int(math.floor(math.log10(microstructure_cropped.max())+1)) + format = '%g' if options.real else '%{}i'.format(int(math.floor(math.log10(microstructure_cropped.max())+1))) table.data = microstructure_cropped.reshape((newInfo['grid'][0],newInfo['grid'][1]*newInfo['grid'][2]),order='F').transpose() - table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ') + table.data_writeArray(format,delimiter=' ') # --- output finalization -------------------------------------------------------------------------- From e39441f36bbff3d2dc4d63e42b2492bb6081c16a Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 17 Mar 2016 19:02:30 -0400 Subject: [PATCH 15/18] corrected usage message to tell about string comparison in filter condition --- processing/post/filterTable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/filterTable.py b/processing/post/filterTable.py index d375785dc..6723d9faf 100755 --- a/processing/post/filterTable.py +++ b/processing/post/filterTable.py @@ -35,7 +35,7 @@ Filter rows according to condition and columns by either white or black listing. Examples: Every odd row if x coordinate is positive -- " #ip.x# >= 0.0 and #_row_#%2 == 1 ). -All rows where label 'foo' equals 'bar' -- " #foo# == \"bar\" " +All rows where label 'foo' equals 'bar' -- " #s#foo# == 'bar' " """, version = scriptID) From 2506667326235ed45fc461160fca211a80e4ef2c Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 18 Mar 2016 04:13:24 +0100 Subject: [PATCH 16/18] updated version information after successful test of v2.0.0-43-ge39441f --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 09b18f22c..690c6cb1d 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-36-g2b7524e +v2.0.0-43-ge39441f From 6865987ea743fc5a2c0f92c42e37766b1ef90dea Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 20 Mar 2016 15:01:47 -0400 Subject: [PATCH 17/18] restored option to output microstructure indices as one-dimensional list (not clear why that got lost!) --- processing/pre/geom_unpack.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/processing/pre/geom_unpack.py b/processing/pre/geom_unpack.py index a07166cdd..85c60d6e5 100755 --- a/processing/pre/geom_unpack.py +++ b/processing/pre/geom_unpack.py @@ -18,6 +18,14 @@ Unpack geometry files containing ranges "a to b" and/or "n of x" multiples (excl """, version = scriptID) +parser.add_option('-1', '--onedimensional', + dest = 'oneD', + action = 'store_true', + help = 'output geom file with one-dimensional data arrangement') + +parser.set_defaults(oneD = False, + ) + (options, filenames) = parser.parse_args() # --- loop over input files ------------------------------------------------------------------------- @@ -69,7 +77,8 @@ for name in filenames: microstructure = table.microstructure_read(info['grid']) # read microstructure formatwidth = int(math.floor(math.log10(microstructure.max())+1)) # efficient number printing format - table.data = microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose() + table.data = microstructure if options.oneD else \ + microstructure.reshape((info['grid'][0],info['grid'][1]*info['grid'][2]),order='F').transpose() table.data_writeArray('%%%ii'%(formatwidth),delimiter = ' ') #--- output finalization -------------------------------------------------------------------------- From a41cd3df4068af0cd752d72cc060f4d91204f470 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 20 Mar 2016 16:46:35 -0400 Subject: [PATCH 18/18] added "strict" option to reading microstructure from geom file. will now complain if data mismatches with grid. --- lib/damask/asciitable.py | 8 +++++--- processing/pre/geom_check.py | 27 +++++++++++++++------------ 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index 5212ef00c..62d38e4dd 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -518,7 +518,8 @@ class ASCIItable(): # ------------------------------------------------------------------ def microstructure_read(self, grid, - type = 'i'): + type = 'i', + strict = False): """read microstructure data (from .geom format)""" def datatype(item): return int(item) if type.lower() == 'i' else float(item) @@ -537,6 +538,7 @@ class ASCIItable(): s = min(len(items), N-i) # prevent overflow of microstructure array microstructure[i:i+s] = items[:s] - i += s + i += len(items) - return microstructure + return (microstructure, i == N and not self.data_read() if strict # check for proper point count and end of file + else microstructure) diff --git a/processing/pre/geom_check.py b/processing/pre/geom_check.py index d607a1f73..518e8c3b0 100755 --- a/processing/pre/geom_check.py +++ b/processing/pre/geom_check.py @@ -54,12 +54,25 @@ for name in filenames: errors = [] if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.') + +#--- read microstructure information -------------------------------------------------------------- + + if options.data: + microstructure,ok = table.microstructure_read(info['grid'],strict = True) # read microstructure + + if ok: + structure = vtk.vtkIntArray() + structure.SetName('Microstructures') + for idx in microstructure: structure.InsertNextValue(idx) + + else: errors.append('mismatch between data and grid dimension.') + if errors != []: damask.util.croak(errors) table.close(dismiss = True) continue -# --- generate VTK rectilinear grid -------------------------------------------------------------------------------- +# --- generate VTK rectilinear grid --------------------------------------------------------------- grid = vtk.vtkRectilinearGrid() grid.SetDimensions([x+1 for x in info['grid']]) @@ -72,18 +85,8 @@ for name in filenames: elif i == 1: grid.SetYCoordinates(temp) elif i == 2: grid.SetZCoordinates(temp) -#--- read microstructure information -------------------------------------------------------------- - if options.data: - microstructure = table.microstructure_read(info['grid']) # read microstructure - - structure = vtk.vtkIntArray() - structure.SetName('Microstructures') - - for idx in microstructure: - structure.InsertNextValue(idx) - - grid.GetCellData().AddArray(structure) + if options.data: grid.GetCellData().AddArray(structure) # --- write data ----------------------------------------------------------------------------------- if name: