Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Martin Diehl 2016-06-29 12:39:14 +02:00
commit fbf67de37a
6 changed files with 42 additions and 44 deletions

View File

@ -1 +1 @@
v2.0.0-292-g2ebc5ec v2.0.0-297-ga27aba1

View File

@ -74,32 +74,30 @@ for name in filenames:
} }
# ------------------------------------------ Evaluate condition --------------------------------------- # ------------------------------------------ Evaluate condition ---------------------------------------
if options.condition: if options.condition is not None:
interpolator = [] interpolator = []
condition = options.condition # copy per file, since might be altered inline condition = options.condition # copy per file, since might be altered inline
breaker = False breaker = False
for position,operand in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups for position,operand in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups
condition = condition.replace('#'+operand[0]+'#', condition = condition.replace('#'+operand[0]+'#',
{ '': '{%i}'%position, { '': '{%i}'%position,
's#':'"{%i}"'%position}[operand[1]]) 's#':'"{%i}"'%position}[operand[1]])
if operand[2] in specials: # special label if operand[2] in specials: # special label
interpolator += ['specials["%s"]'%operand[2]] interpolator += ['specials["%s"]'%operand[2]]
else: else:
try: try:
interpolator += ['%s(table.data[%i])'%({ '':'float', interpolator += ['%s(table.data[%i])'%({ '':'float',
's#':'str'}[operand[1]], 's#':'str'}[operand[1]],
table.label_index(operand[2]))] # ccould be generalized to indexrange as array lookup table.label_index(operand[2]))] # could be generalized to indexrange as array lookup
except: except:
damask.util.croak('column "{}" not found.'.format(operand[2])) damask.util.croak('column "{}" not found.'.format(operand[2]))
breaker = True breaker = True
if breaker: continue # found mistake in condition evaluation --> next file if breaker: continue # found mistake in condition evaluation --> next file
evaluator_condition = "'" + condition + "'.format(" + ','.join(interpolator) + ")" evaluator_condition = "'" + condition + "'.format(" + ','.join(interpolator) + ")"
else: condition = ''
# ------------------------------------------ build formulae ---------------------------------------- # ------------------------------------------ build formulae ----------------------------------------
evaluator = {} evaluator = {}
@ -165,19 +163,19 @@ for name in filenames:
for label in output.labels(): for label in output.labels():
oldIndices = table.label_indexrange(label) oldIndices = table.label_indexrange(label)
Nold = max(1,len(oldIndices)) # Nold could be zero for new columns Nold = max(1,len(oldIndices)) # Nold could be zero for new columns
Nnew = len(output.label_indexrange(label)) Nnew = len(output.label_indexrange(label))
output.data_append(eval(evaluator[label]) if label in options.labels and output.data_append(eval(evaluator[label]) if label in options.labels and
(condition == '' or eval(eval(evaluator_condition))) (options.condition is None or eval(eval(evaluator_condition)))
else np.tile([table.data[i] for i in oldIndices] else np.tile([table.data[i] for i in oldIndices]
if label in tabLabels if label in tabLabels
else np.nan, else np.nan,
np.ceil(float(Nnew)/Nold))[:Nnew]) # spread formula result into given number of columns np.ceil(float(Nnew)/Nold))[:Nnew]) # spread formula result into given number of columns
outputAlive = output.data_write() # output processed line outputAlive = output.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------
table.input_close() # close ASCII tables table.input_close() # close ASCII tables
output.close() # close ASCII tables output.close() # close ASCII tables

View File

@ -18,7 +18,7 @@ def curlFFT(geomdim,field):
if n == 3: dataType = 'vector' if n == 3: dataType = 'vector'
elif n == 9: dataType = 'tensor' elif n == 9: dataType = 'tensor'
field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
curl_fourier = np.empty(field_fourier.shape,'c16') curl_fourier = np.empty(field_fourier.shape,'c16')
# differentiation in Fourier space # differentiation in Fourier space
@ -56,7 +56,7 @@ def curlFFT(geomdim,field):
curl_fourier[i,j,k,2] = ( field_fourier[i,j,k,1]*xi[0]\ curl_fourier[i,j,k,2] = ( field_fourier[i,j,k,1]*xi[0]\
-field_fourier[i,j,k,0]*xi[1]) *TWOPIIMG -field_fourier[i,j,k,0]*xi[1]) *TWOPIIMG
return np.fft.fftpack.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n]) return np.fft.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -158,7 +158,7 @@ for name in filenames:
# 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], stack.append(curlFFT(size[::-1],
table.data[:,data['column'][i]:data['column'][i]+data['dim']]. table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape([grid[2],grid[1],grid[0]]+data['shape']))) reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result ----------------------------------------- # ------------------------------------------ output result -----------------------------------------

View File

@ -18,7 +18,7 @@ def cell2node(cellData,grid):
datalen = np.array(cellData.shape[3:]).prod() datalen = np.array(cellData.shape[3:]).prod()
for i in xrange(datalen): for i in xrange(datalen):
node = scipy.ndimage.convolve(cellData.reshape(tuple(grid)+(datalen,))[...,i], node = scipy.ndimage.convolve(cellData.reshape(tuple(grid[::-1])+(datalen,))[...,i],
np.ones((2,2,2))/8., # 2x2x2 neighborhood of cells np.ones((2,2,2))/8., # 2x2x2 neighborhood of cells
mode = 'wrap', mode = 'wrap',
origin = -1, # offset to have cell origin as center origin = -1, # offset to have cell origin as center
@ -35,14 +35,14 @@ def cell2node(cellData,grid):
def displacementAvgFFT(F,grid,size,nodal=False,transformed=False): 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""" """calculate average cell center (or nodal) displacement for deformation gradient field specified in each grid cell"""
if nodal: if nodal:
x, y, z = np.meshgrid(np.linspace(0,size[0],1+grid[0]), 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[1],1+grid[1]),
np.linspace(0,size[2],1+grid[2]), np.linspace(0,size[0],1+grid[0]),
indexing = 'ij') indexing = 'ij')
else: else:
x, y, z = np.meshgrid(np.linspace(0,size[0],grid[0],endpoint=False), x, y, z = np.meshgrid(np.linspace(0,size[2],grid[2],endpoint=False),
np.linspace(0,size[1],grid[1],endpoint=False), np.linspace(0,size[1],grid[1],endpoint=False),
np.linspace(0,size[2],grid[2],endpoint=False), np.linspace(0,size[0],grid[0],endpoint=False),
indexing = 'ij') indexing = 'ij')
origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
@ -69,7 +69,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# integration in Fourier space # integration in Fourier space
displacement_fourier = +np.einsum('ijkml,ijkl,l->ijkm', displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm',
F if transformed else np.fft.rfftn(F,axes=(0,1,2)), F if transformed else np.fft.rfftn(F,axes=(0,1,2)),
k_s, k_s,
integrator, integrator,
@ -78,7 +78,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
#-------------------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------------------
# backtransformation to real space # backtransformation to real space
displacement = np.fft.irfftn(displacement_fourier,grid,axes=(0,1,2)) displacement = np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2))
return cell2node(displacement,grid) if nodal else displacement return cell2node(displacement,grid) if nodal else displacement
@ -186,8 +186,8 @@ for name in filenames:
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... 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...
displacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True) fluctDisplacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True)
avgDisplacement = displacementAvgFFT (F_fourier,grid,size,options.nodal,transformed=True) avgDisplacement = displacementAvgFFT (F_fourier,grid,size,options.nodal,transformed=True)
# ------------------------------------------ assemble header --------------------------------------- # ------------------------------------------ assemble header ---------------------------------------
@ -203,18 +203,18 @@ for name in filenames:
# ------------------------------------------ output data ------------------------------------------- # ------------------------------------------ output data -------------------------------------------
zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else xrange(grid[2]) Zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else xrange(grid[2])
yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else xrange(grid[1]) Yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else xrange(grid[1])
xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0]) Xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0])
for i,z in enumerate(zrange): for i,z in enumerate(Zrange):
for j,y in enumerate(yrange): for j,y in enumerate(Yrange):
for k,x in enumerate(xrange): for k,x in enumerate(Xrange):
if options.nodal: table.data_clear() if options.nodal: table.data_clear()
else: table.data_read() else: table.data_read()
table.data_append([x,y,z] if options.nodal else []) table.data_append([x,y,z] if options.nodal else [])
table.data_append(list(avgDisplacement[i,j,k,:])) table.data_append(list( avgDisplacement[i,j,k,:]))
table.data_append(list( displacement[i,j,k,:])) table.data_append(list(fluctDisplacement[i,j,k,:]))
table.data_write() table.data_write()
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

View File

@ -15,7 +15,7 @@ def divFFT(geomdim,field):
N = grid.prod() # field size N = grid.prod() # field size
n = np.array(np.shape(field)[3:]).prod() # data size n = np.array(np.shape(field)[3:]).prod() # data size
field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector div_fourier = np.empty(field_fourier.shape[0:len(np.shape(field))-1],'c16') # size depents on whether tensor or vector
# differentiation in Fourier space # differentiation in Fourier space
@ -42,7 +42,7 @@ def divFFT(geomdim,field):
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 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),s=shapeFFT).reshape([N,n/3]) return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n/3])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -145,7 +145,7 @@ for name in filenames:
# 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[:,data['column'][i]:data['column'][i]+data['dim']]. table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape([grid[2],grid[1],grid[0]]+data['shape']))) reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result ----------------------------------------- # ------------------------------------------ output result -----------------------------------------

View File

@ -18,7 +18,7 @@ def gradFFT(geomdim,field):
if n == 3: dataType = 'vector' if n == 3: dataType = 'vector'
elif n == 1: dataType = 'scalar' elif n == 1: dataType = 'scalar'
field_fourier = np.fft.fftpack.rfftn(field,axes=(0,1,2),s=shapeFFT) field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT)
grad_fourier = np.empty(field_fourier.shape+(3,),'c16') grad_fourier = np.empty(field_fourier.shape+(3,),'c16')
# differentiation in Fourier space # differentiation in Fourier space
@ -46,7 +46,7 @@ def gradFFT(geomdim,field):
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 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),s=shapeFFT).reshape([N,3*n]) return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n])
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -148,7 +148,7 @@ for name in filenames:
# 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], stack.append(gradFFT(size[::-1],
table.data[:,data['column'][i]:data['column'][i]+data['dim']]. table.data[:,data['column'][i]:data['column'][i]+data['dim']].
reshape([grid[2],grid[1],grid[0]]+data['shape']))) reshape(grid[::-1].tolist()+data['shape'])))
# ------------------------------------------ output result ----------------------------------------- # ------------------------------------------ output result -----------------------------------------