functions for spatial coordinates on regular grids

This commit is contained in:
Martin Diehl 2019-12-03 18:59:59 +01:00
parent 62ca2952fc
commit e006e0ebec
1 changed files with 45 additions and 22 deletions

View File

@ -2,7 +2,7 @@ import numpy as np
def __ks(size,field,first_order=False): def __ks(size,field,first_order=False):
"""Get wave numbers operator.""" """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] 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) 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 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 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): 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 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 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): 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 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 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): def coord0_cell(grid,size):
"""Positions of nodes (undeformed).""" """Cell center 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 coord_cell(grid,size):
"""Positions of cell centers (undeformed)."""
delta = size/grid*0.5 delta = size/grid*0.5
x, y, z = np.meshgrid(np.linspace(delta[2],size[2]-delta[2],grid[2]), 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]), 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) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def displacement_fluct_cell(size,F):
def displacement_fluct(size,F): """Cell center displacement field from fluctuation part of the deformation gradient field."""
"""Calculate displacement field from deformation gradient field.""" integrator = 0.5j*size/np.pi
integrator = 0.5j * size / np.pi
k_s = __ks(size,F,False) 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', displacement = -np.einsum('ijkml,ijkl,l->ijkm',
np.fft.rfftn(F,axes=(0,1,2)), np.fft.rfftn(F,axes=(0,1,2)),
k_s, k_s,
integrator, 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')