diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index e49ff47a9..bf757bdb9 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -2,7 +2,7 @@ import numpy as np def __ks(size,field,first_order=False): """Get wave numbers operator.""" - grid = np.array(np.shape(field)[0:3]) + grid = np.array(np.shape(field)[:3]) k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0] if grid[0]%2 == 0 and first_order: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) @@ -29,7 +29,7 @@ def curl(size,field): curl = (np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 3 np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3x3 - return np.fft.irfftn(curl,axes=(0,1,2),s=field.shape[0:3]) + return np.fft.irfftn(curl,axes=(0,1,2),s=field.shape[:3]) def divergence(size,field): @@ -41,7 +41,7 @@ def divergence(size,field): divergence = (np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*np.pi if n == 3 else # vector, 3 -> 1 np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 3 - return np.fft.irfftn(divergence,axes=(0,1,2),s=field.shape[0:3]) + return np.fft.irfftn(divergence,axes=(0,1,2),s=field.shape[:3]) def gradient(size,field): @@ -53,21 +53,11 @@ def gradient(size,field): gradient = (np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*np.pi if n == 1 else # scalar, 1 -> 3 np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*np.pi) # vector, 3 -> 3x3 - return np.fft.irfftn(gradient,axes=(0,1,2),s=field.shape[0:3]) + return np.fft.irfftn(gradient,axes=(0,1,2),s=field.shape[:3]) -def coord_node(grid,size): - """Positions of nodes (undeformed).""" - 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') - - return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) - - -def coord_cell(grid,size): - """Positions of cell centers (undeformed).""" +def coord0_cell(grid,size): + """Cell center positions (undeformed).""" 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]), @@ -76,17 +66,50 @@ def coord_cell(grid,size): return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) - -def displacement_fluct(size,F): - """Calculate displacement field from deformation gradient field.""" - integrator = 0.5j * size / np.pi +def displacement_fluct_cell(size,F): + """Cell center displacement field from fluctuation part of the deformation gradient field.""" + integrator = 0.5j*size/np.pi k_s = __ks(size,F,False) + k_s_squared = np.einsum('...l,...l',k_s,k_s) + k_s_squared[0,0,0] = 1.0 displacement = -np.einsum('ijkml,ijkl,l->ijkm', np.fft.rfftn(F,axes=(0,1,2)), k_s, integrator, - ) / k_sSquared[...,np.newaxis] + ) / k_s_squared[...,np.newaxis] - return np.fft.irfftn(displacement,axes=(0,1,2)) + return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) + +def displacement_avg_cell(size,F): + """Cell center displacement field from average part of the deformation gradient field.""" + F_avg = np.average(F,axis=(0,1,2)) + return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),coord0_cell(F.shape[:3],size)) + + +def coord0_node(grid,size): + """Nodal positions (undeformed).""" + 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') + + return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) + +def displacement_fluct_node(size,F): + return cell_2_node(displacement_fluct_cell(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)) + + +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 + + return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')