using 3 way merge to have syntax as similar as possible
This commit is contained in:
parent
cbe2fb8d57
commit
873c52cceb
|
@ -12,8 +12,8 @@ scriptID = ' '.join([scriptName,damask.version])
|
||||||
def curlFFT(geomdim,field):
|
def curlFFT(geomdim,field):
|
||||||
shapeFFT = np.array(np.shape(field))[0:3]
|
shapeFFT = np.array(np.shape(field))[0:3]
|
||||||
grid = np.array(np.shape(field)[2::-1])
|
grid = np.array(np.shape(field)[2::-1])
|
||||||
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
|
||||||
|
|
||||||
if n == 3: dataType = 'vector'
|
if n == 3: dataType = 'vector'
|
||||||
elif n == 9: dataType = 'tensor'
|
elif n == 9: dataType = 'tensor'
|
||||||
|
@ -24,10 +24,10 @@ def curlFFT(geomdim,field):
|
||||||
# differentiation in Fourier space
|
# differentiation in Fourier space
|
||||||
TWOPIIMG = 2.0j*math.pi
|
TWOPIIMG = 2.0j*math.pi
|
||||||
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
||||||
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||||
|
|
||||||
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
||||||
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||||
|
|
||||||
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
||||||
|
|
||||||
|
@ -35,12 +35,12 @@ def curlFFT(geomdim,field):
|
||||||
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
||||||
|
|
||||||
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 symbols
|
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = 1.0 # Levi-Civita symbols
|
||||||
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
|
||||||
|
|
||||||
if dataType == 'tensor': # tensor, 3x3 -> 3x3
|
if dataType == 'tensor': # tensor, 3x3 -> 3x3
|
||||||
curl_fourier = np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*TWOPIIMG
|
curl_fourier = np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*TWOPIIMG
|
||||||
elif dataType == 'vector': # vector, 3 -> 3
|
elif dataType == 'vector': # vector, 3 -> 3
|
||||||
curl_fourier = np.einsum('slm,ijkl,ijkm->ijks',e,k_s,field_fourier)*TWOPIIMG
|
curl_fourier = np.einsum('slm,ijkl,ijkm->ijks',e,k_s,field_fourier)*TWOPIIMG
|
||||||
|
|
||||||
return np.fft.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])
|
||||||
|
|
|
@ -12,30 +12,30 @@ scriptID = ' '.join([scriptName,damask.version])
|
||||||
def divFFT(geomdim,field):
|
def divFFT(geomdim,field):
|
||||||
shapeFFT = np.array(np.shape(field))[0:3]
|
shapeFFT = np.array(np.shape(field))[0:3]
|
||||||
grid = np.array(np.shape(field)[2::-1])
|
grid = np.array(np.shape(field)[2::-1])
|
||||||
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
|
||||||
|
|
||||||
if n == 3: dataType = 'vector'
|
if n == 3: dataType = 'vector'
|
||||||
elif n == 9: dataType = 'tensor'
|
elif n == 9: dataType = 'tensor'
|
||||||
|
|
||||||
field_fourier = np.fft.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')
|
||||||
|
|
||||||
# differentiation in Fourier space
|
# differentiation in Fourier space
|
||||||
TWOPIIMG = 2.0j*math.pi
|
TWOPIIMG = 2.0j*math.pi
|
||||||
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
||||||
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||||
|
|
||||||
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
||||||
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||||
|
|
||||||
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
||||||
|
|
||||||
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
||||||
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
||||||
if dataType == 'tensor': # tensor, 3x3 -> 3
|
if dataType == 'tensor': # tensor, 3x3 -> 3
|
||||||
div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG
|
div_fourier = np.einsum('ijklm,ijkm->ijkl',field_fourier,k_s)*TWOPIIMG
|
||||||
elif dataType == 'vector': # vector, 3 -> 1
|
elif dataType == 'vector': # vector, 3 -> 1
|
||||||
div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG
|
div_fourier = np.einsum('ijkl,ijkl->ijk',field_fourier,k_s)*TWOPIIMG
|
||||||
|
|
||||||
return np.fft.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])
|
||||||
|
|
|
@ -12,8 +12,8 @@ scriptID = ' '.join([scriptName,damask.version])
|
||||||
def gradFFT(geomdim,field):
|
def gradFFT(geomdim,field):
|
||||||
shapeFFT = np.array(np.shape(field))[0:3]
|
shapeFFT = np.array(np.shape(field))[0:3]
|
||||||
grid = np.array(np.shape(field)[2::-1])
|
grid = np.array(np.shape(field)[2::-1])
|
||||||
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
|
||||||
|
|
||||||
if n == 3: dataType = 'vector'
|
if n == 3: dataType = 'vector'
|
||||||
elif n == 1: dataType = 'scalar'
|
elif n == 1: dataType = 'scalar'
|
||||||
|
@ -24,18 +24,18 @@ def gradFFT(geomdim,field):
|
||||||
# differentiation in Fourier space
|
# differentiation in Fourier space
|
||||||
TWOPIIMG = 2.0j*math.pi
|
TWOPIIMG = 2.0j*math.pi
|
||||||
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/geomdim[0]
|
||||||
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||||
|
|
||||||
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/geomdim[1]
|
||||||
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # for even grid, set Nyquist freq to 0 (Johnson, MIT, 2011)
|
||||||
|
|
||||||
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
k_si = np.arange(grid[0]//2+1)/geomdim[2]
|
||||||
|
|
||||||
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij')
|
||||||
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')
|
||||||
if dataType == 'vector': # vector, 3 -> 3x3
|
if dataType == 'vector': # vector, 3 -> 3x3
|
||||||
grad_fourier = np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*TWOPIIMG
|
grad_fourier = np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*TWOPIIMG
|
||||||
elif dataType == 'scalar': # scalar, 1 -> 3
|
elif dataType == 'scalar': # scalar, 1 -> 3
|
||||||
grad_fourier = np.einsum('ijkl,ijkl->ijkl',field_fourier,k_s)*TWOPIIMG
|
grad_fourier = np.einsum('ijkl,ijkl->ijkl',field_fourier,k_s)*TWOPIIMG
|
||||||
|
|
||||||
return np.fft.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])
|
||||||
|
|
Loading…
Reference in New Issue