proper spacing

This commit is contained in:
Martin Diehl 2020-02-20 22:50:54 +01:00
parent ab88ffd28e
commit a670ab269f
1 changed files with 41 additions and 29 deletions

View File

@ -1,14 +1,14 @@
from scipy import spatial from scipy import spatial
import numpy as np import numpy as np
def __ks(size,grid,first_order=False): def _ks(size,grid,first_order=False):
""" """
Get wave numbers operator. Get wave numbers operator.
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
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]
@ -30,14 +30,14 @@ def curl(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
e = np.zeros((3, 3, 3)) e = np.zeros((3, 3, 3))
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol
e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0 e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
@ -54,11 +54,11 @@ def divergence(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
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
@ -74,11 +74,11 @@ def gradient(size,field):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
""" """
n = np.prod(field.shape[3:]) n = np.prod(field.shape[3:])
k_s = __ks(size,field.shape[:3],True) k_s = _ks(size,field.shape[:3],True)
field_fourier = np.fft.rfftn(field,axes=(0,1,2)) field_fourier = np.fft.rfftn(field,axes=(0,1,2))
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
@ -96,7 +96,7 @@ def cell_coord0(grid,size,origin=np.zeros(3)):
grid : numpy.ndarray grid : numpy.ndarray
number of grid points. number of grid points.
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
physical origin of the periodic field. Default is [0.0,0.0,0.0]. physical origin of the periodic field. Default is [0.0,0.0,0.0].
@ -108,7 +108,8 @@ def cell_coord0(grid,size,origin=np.zeros(3)):
np.linspace(start[0],end[0],grid[0]), np.linspace(start[0],end[0],grid[0]),
indexing = 'ij') indexing = 'ij')
return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def cell_displacement_fluct(size,F): def cell_displacement_fluct(size,F):
""" """
@ -117,14 +118,14 @@ def cell_displacement_fluct(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
integrator = 0.5j*size/np.pi integrator = 0.5j*size/np.pi
k_s = __ks(size,F.shape[:3],False) k_s = _ks(size,F.shape[:3],False)
k_s_squared = np.einsum('...l,...l',k_s,k_s) k_s_squared = np.einsum('...l,...l',k_s,k_s)
k_s_squared[0,0,0] = 1.0 k_s_squared[0,0,0] = 1.0
@ -136,6 +137,7 @@ def cell_displacement_fluct(size,F):
return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3])
def cell_displacement_avg(size,F): def cell_displacement_avg(size,F):
""" """
Cell center displacement field from average part of the deformation gradient field. Cell center displacement field from average part of the deformation gradient field.
@ -143,7 +145,7 @@ def cell_displacement_avg(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -151,6 +153,7 @@ def cell_displacement_avg(size,F):
F_avg = np.average(F,axis=(0,1,2)) F_avg = np.average(F,axis=(0,1,2))
return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3][::-1],size)) return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3][::-1],size))
def cell_displacement(size,F): def cell_displacement(size,F):
""" """
Cell center displacement field from deformation gradient field. Cell center displacement field from deformation gradient field.
@ -158,13 +161,14 @@ def cell_displacement(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F) return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F)
def cell_coord(size,F,origin=np.zeros(3)): def cell_coord(size,F,origin=np.zeros(3)):
""" """
Cell center positions. Cell center positions.
@ -172,7 +176,7 @@ def cell_coord(size,F,origin=np.zeros(3)):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
@ -181,6 +185,7 @@ def cell_coord(size,F,origin=np.zeros(3)):
""" """
return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F) return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F)
def cell_coord0_gridSizeOrigin(coord0,ordered=True): def cell_coord0_gridSizeOrigin(coord0,ordered=True):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of cell positions. Return grid 'DNA', i.e. grid, size, and origin from array of cell positions.
@ -200,11 +205,11 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
size = grid/np.maximum(grid-1,1) * (maxcorner-mincorner) size = grid/np.maximum(grid-1,1) * (maxcorner-mincorner)
delta = size/grid delta = size/grid
origin = mincorner - delta*.5 origin = mincorner - delta*.5
# 1D/2D: size/origin combination undefined, set origin to 0.0 # 1D/2D: size/origin combination undefined, set origin to 0.0
size [np.where(grid==1)] = origin[np.where(grid==1)]*2. size [np.where(grid==1)] = origin[np.where(grid==1)]*2.
origin[np.where(grid==1)] = 0.0 origin[np.where(grid==1)] = 0.0
if grid.prod() != len(coord0): if grid.prod() != len(coord0):
raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))
@ -221,6 +226,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True):
return (grid,size,origin) return (grid,size,origin)
def coord0_check(coord0): def coord0_check(coord0):
""" """
Check whether coordinates lie on a regular grid. Check whether coordinates lie on a regular grid.
@ -234,7 +240,6 @@ def coord0_check(coord0):
cell_coord0_gridSizeOrigin(coord0,ordered=True) cell_coord0_gridSizeOrigin(coord0,ordered=True)
def node_coord0(grid,size,origin=np.zeros(3)): def node_coord0(grid,size,origin=np.zeros(3)):
""" """
Nodal positions (undeformed). Nodal positions (undeformed).
@ -244,7 +249,7 @@ def node_coord0(grid,size,origin=np.zeros(3)):
grid : numpy.ndarray grid : numpy.ndarray
number of grid points. number of grid points.
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
physical origin of the periodic field. Default is [0.0,0.0,0.0]. physical origin of the periodic field. Default is [0.0,0.0,0.0].
@ -253,8 +258,9 @@ def node_coord0(grid,size,origin=np.zeros(3)):
np.linspace(origin[1],size[1]+origin[1],1+grid[1]), np.linspace(origin[1],size[1]+origin[1],1+grid[1]),
np.linspace(origin[0],size[0]+origin[0],1+grid[0]), np.linspace(origin[0],size[0]+origin[0],1+grid[0]),
indexing = 'ij') indexing = 'ij')
return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3)
def node_displacement_fluct(size,F): def node_displacement_fluct(size,F):
""" """
@ -263,13 +269,14 @@ def node_displacement_fluct(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
return cell_2_node(cell_displacement_fluct(size,F)) return cell_2_node(cell_displacement_fluct(size,F))
def node_displacement_avg(size,F): def node_displacement_avg(size,F):
""" """
Nodal displacement field from average part of the deformation gradient field. Nodal displacement field from average part of the deformation gradient field.
@ -277,7 +284,7 @@ def node_displacement_avg(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
@ -285,6 +292,7 @@ def node_displacement_avg(size,F):
F_avg = np.average(F,axis=(0,1,2)) F_avg = np.average(F,axis=(0,1,2))
return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3][::-1],size)) return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3][::-1],size))
def node_displacement(size,F): def node_displacement(size,F):
""" """
Nodal displacement field from deformation gradient field. Nodal displacement field from deformation gradient field.
@ -292,13 +300,14 @@ def node_displacement(size,F):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
""" """
return node_displacement_avg(size,F) + node_displacement_fluct(size,F) return node_displacement_avg(size,F) + node_displacement_fluct(size,F)
def node_coord(size,F,origin=np.zeros(3)): def node_coord(size,F,origin=np.zeros(3)):
""" """
Nodal positions. Nodal positions.
@ -306,7 +315,7 @@ def node_coord(size,F,origin=np.zeros(3)):
Parameters Parameters
---------- ----------
size : numpy.ndarray size : numpy.ndarray
physical size of the periodic field. physical size of the periodic field.
F : numpy.ndarray F : numpy.ndarray
deformation gradient field. deformation gradient field.
origin : numpy.ndarray, optional origin : numpy.ndarray, optional
@ -315,22 +324,25 @@ def node_coord(size,F,origin=np.zeros(3)):
""" """
return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F) return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F)
def cell_2_node(cell_data): def cell_2_node(cell_data):
"""Interpolate periodic cell data to nodal data.""" """Interpolate periodic cell data to nodal data."""
n = ( cell_data + np.roll(cell_data,1,(0,1,2)) 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,)) + 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 + 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') return np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_2_cell(node_data): def node_2_cell(node_data):
"""Interpolate periodic nodal data to cell data.""" """Interpolate periodic nodal data to cell data."""
c = ( node_data + np.roll(node_data,1,(0,1,2)) c = ( node_data + np.roll(node_data,1,(0,1,2))
+ np.roll(node_data,1,(0,)) + np.roll(node_data,1,(1,)) + np.roll(node_data,1,(2,)) + np.roll(node_data,1,(0,)) + np.roll(node_data,1,(1,)) + np.roll(node_data,1,(2,))
+ np.roll(node_data,1,(0,1)) + np.roll(node_data,1,(1,2)) + np.roll(node_data,1,(2,0)))*0.125 + np.roll(node_data,1,(0,1)) + np.roll(node_data,1,(1,2)) + np.roll(node_data,1,(2,0)))*0.125
return c[:-1,:-1,:-1] return c[:-1,:-1,:-1]
def node_coord0_gridSizeOrigin(coord0,ordered=False): def node_coord0_gridSizeOrigin(coord0,ordered=False):
""" """
Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions. Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions.
@ -349,7 +361,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False):
grid = np.array(list(map(len,coords)),'i') - 1 grid = np.array(list(map(len,coords)),'i') - 1
size = maxcorner-mincorner size = maxcorner-mincorner
origin = mincorner origin = mincorner
if (grid+1).prod() != len(coord0): if (grid+1).prod() != len(coord0):
raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid))