fixed fluct displacement

This commit is contained in:
Martin Diehl 2016-06-29 12:11:07 +02:00
parent a64388355a
commit 9bf8b6221c
1 changed files with 12 additions and 12 deletions

View File

@ -18,7 +18,7 @@ def cell2node(cellData,grid):
datalen = np.array(cellData.shape[3:]).prod()
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
mode = 'wrap',
origin = -1, # offset to have cell origin as center
@ -69,7 +69,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
#--------------------------------------------------------------------------------------------------
# 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)),
k_s,
integrator,
@ -78,7 +78,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
#--------------------------------------------------------------------------------------------------
# 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
@ -187,7 +187,7 @@ 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...
displacement = 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 ---------------------------------------
@ -196,24 +196,24 @@ for name in filenames:
table.labels_clear()
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
#['{}_avg({}).{}' .format(i+1,options.defgrad,options.pos) for i in xrange(3)] +
table.labels_append((['{}_pos' .format(i+1) for i in xrange(3)] if options.nodal else []) +
['{}_avg({}).{}' .format(i+1,options.defgrad,options.pos) for i in xrange(3)] +
['{}_fluct({}).{}'.format(i+1,options.defgrad,options.pos) for i in xrange(3)] )
table.head_write()
# ------------------------------------------ output data -------------------------------------------
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])
xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0])
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])
Xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else xrange(grid[0])
for i,z in enumerate(zrange):
for j,y in enumerate(yrange):
for k,x in enumerate(xrange):
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(avgDisplacement[i,j,k,:]))
table.data_append(list( displacement[i,j,k,:]))
table.data_write()