From 23f5e0fa58e8ca0137236bd8499a53534f11a2cb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 26 Nov 2019 10:25:39 +0100 Subject: [PATCH 01/36] filters for operations on regular grids (in fourier space) --- python/damask/grid_filters.py | 113 ++++++++++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) create mode 100644 python/damask/grid_filters.py diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py new file mode 100644 index 000000000..6588e59bc --- /dev/null +++ b/python/damask/grid_filters.py @@ -0,0 +1,113 @@ +import numpy as np + + +def curl(size,field): + """Calculate curl of a vector or tensor field in Fourier space.""" + shapeFFT = np.array(np.shape(field))[0:3] + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size + + field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) + curl_fourier = np.empty(field_fourier.shape,'c16') + + k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] + if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + + k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] + if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + + k_si = np.arange(grid[0]//2+1)/size[2] + + 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') + + 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, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0 + + curl_fourier = 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_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n]) + + +def divergence(size,field): + """Calculate divergence of a vector or tensor field in Fourier space.""" + shapeFFT = np.array(np.shape(field))[0:3] + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size + + 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') + + k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] + if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + + k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] + if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + + k_si = np.arange(grid[0]//2+1)/size[2] + + 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') + + div_fourier = 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(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n//3]) + + +def gradient(size,field): + """Calculate gradient of a vector or scalar field in Fourier space.""" + shapeFFT = np.array(np.shape(field))[0:3] + grid = np.array(np.shape(field)[2::-1]) + N = grid.prod() # field size + n = np.array(np.shape(field)[3:]).prod() # data size + + field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) + grad_fourier = np.empty(field_fourier.shape+(3,),'c16') + + k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] + if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + + k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] + if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + + k_si = np.arange(grid[0]//2+1)/size[2] + + 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') + grad_fourier = 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(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n]) + + +#-------------------------------------------------------------------------------------------------- +def displacementFluctFFT(F,size): + """Calculate displacement field from deformation gradient field.""" + integrator = 0.5j * size / np.pi + + kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])), + np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])), + np.arange(grid[0]//2+1), + indexing = 'ij') + k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) + k_sSquared = np.einsum('...l,...l',k_s,k_s) + k_sSquared[0,0,0] = 1.0 # ignore global average frequency + +#-------------------------------------------------------------------------------------------------- +# integration in Fourier space + + displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm', + np.fft.rfftn(F,axes=(0,1,2)), + k_s, + integrator, + ) / k_sSquared[...,np.newaxis] + +#-------------------------------------------------------------------------------------------------- +# backtransformation to real space + + return np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2)) From b85049cb811b70d71d497001af07b7e73b4fbc25 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Nov 2019 05:41:53 +0100 Subject: [PATCH 02/36] use brackets for line continuation with comments --- python/damask/__init__.py | 1 + python/damask/grid_filters.py | 15 ++++++++------- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 2cd26cf1c..77666561c 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -23,4 +23,5 @@ from .util import extendableOption # noqa # functions in modules from . import mechanics # noqa +from . import grid_filters # noqa diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 6588e59bc..26ee2cfcc 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -26,10 +26,10 @@ def curl(size,field): 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 - curl_fourier = 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 + 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_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n]) + return np.fft.irfftn(curl,axes=(0,1,2),s=shapeFFT).reshape([N,n]) def divergence(size,field): @@ -53,8 +53,8 @@ def divergence(size,field): 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') - div_fourier = 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 + 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(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n//3]) @@ -79,8 +79,9 @@ def gradient(size,field): 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') - grad_fourier = 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 + + 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(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n]) From 4c4ccfe72e88053822c83886089b96bfda0f5a35 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Nov 2019 06:27:19 +0100 Subject: [PATCH 03/36] not needed --- processing/post/addCurl.py | 1 - processing/post/addDivergence.py | 1 - processing/post/addGradient.py | 1 - 3 files changed, 3 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 484af9677..b4dd465a9 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -27,7 +27,6 @@ def curlFFT(geomdim,field): n = np.array(np.shape(field)[3:]).prod() # data size field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - curl_fourier = np.empty(field_fourier.shape,'c16') # differentiation in Fourier space TWOPIIMG = 2.0j*np.pi diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 31a18f8e1..2d6af2036 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -27,7 +27,6 @@ def divFFT(geomdim,field): n = np.array(np.shape(field)[3:]).prod() # data size 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') # differentiation in Fourier space TWOPIIMG = 2.0j*np.pi diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index bfadb578e..252b72eb8 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -27,7 +27,6 @@ def gradFFT(geomdim,field): n = np.array(np.shape(field)[3:]).prod() # data size field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - grad_fourier = np.empty(field_fourier.shape+(3,),'c16') # differentiation in Fourier space TWOPIIMG = 2.0j*np.pi From 80b50f460e24ad6ef7a49ac2ac88a5d2946aa74e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Nov 2019 10:09:22 +0100 Subject: [PATCH 04/36] cleaning trying to get rid of strange re-ordering related to ASCII table data layout --- python/damask/grid_filters.py | 47 ++++++++++++++--------------------- 1 file changed, 19 insertions(+), 28 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 26ee2cfcc..69ee85033 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -3,21 +3,18 @@ import numpy as np def curl(size,field): """Calculate curl of a vector or tensor field in Fourier space.""" - shapeFFT = np.array(np.shape(field))[0:3] - grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size + grid = np.array(np.shape(field)[0:3]) n = np.array(np.shape(field)[3:]).prod() # data size - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - curl_fourier = np.empty(field_fourier.shape,'c16') + field_fourier = np.fft.rfftn(field,axes=(0,1,2)) - k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] - if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + 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: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - k_si = np.arange(grid[0]//2+1)/size[2] + k_si = np.arange(grid[2]//2+1)/size[2] 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') @@ -26,29 +23,26 @@ def curl(size,field): 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 - 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 + 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=shapeFFT).reshape([N,n]) + return np.fft.irfftn(curl,axes=(0,1,2)) def divergence(size,field): """Calculate divergence of a vector or tensor field in Fourier space.""" - shapeFFT = np.array(np.shape(field))[0:3] - grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size + grid = np.array(np.shape(field)[0:3]) n = np.array(np.shape(field)[3:]).prod() # data size - 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') + field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=grid) - k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] - if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + 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: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - k_si = np.arange(grid[0]//2+1)/size[2] + k_si = np.arange(grid[2]//2+1)/size[2] 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') @@ -56,26 +50,23 @@ 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(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n//3]) + return np.fft.irfftn(div_fourier,axes=(0,1,2),s=grid) def gradient(size,field): """Calculate gradient of a vector or scalar field in Fourier space.""" - shapeFFT = np.array(np.shape(field))[0:3] grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size n = np.array(np.shape(field)[3:]).prod() # data size - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - grad_fourier = np.empty(field_fourier.shape+(3,),'c16') + field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=grid) - k_sk = np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2]))/size[0] - if grid[2]%2 == 0: k_sk[grid[2]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + 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: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - k_si = np.arange(grid[0]//2+1)/size[2] + k_si = np.arange(grid[2]//2+1)/size[2] 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') @@ -83,7 +74,7 @@ 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(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n]) + return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=grid) #-------------------------------------------------------------------------------------------------- From 3e65d44e073f925c4bbd8ed0c53b4a5308c1e8ee Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Nov 2019 15:46:22 +0100 Subject: [PATCH 05/36] centralized facilities for differential operations note the need to reverse the grid shape in data from the ASCII table. If x is fastest, z is slowest we require x to be the rightmost index --- processing/post/addCurl.py | 138 ++++-------------------------- processing/post/addDivergence.py | 140 +++++-------------------------- processing/post/addGradient.py | 135 ++++------------------------- python/damask/grid_filters.py | 63 +++++--------- 4 files changed, 74 insertions(+), 402 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index b4dd465a9..2fcd107c0 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np @@ -12,47 +13,6 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -def merge_dicts(*dict_args): - """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts.""" - result = {} - for dictionary in dict_args: - result.update(dictionary) - return result - -def curlFFT(geomdim,field): - """Calculate curl of a vector or tensor field by transforming into Fourier space.""" - shapeFFT = np.array(np.shape(field))[0:3] - grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size - - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - - # differentiation in Fourier space - TWOPIIMG = 2.0j*np.pi - einsums = { - 3:'slm,ijkl,ijkm->ijks', # vector, 3 -> 3 - 9:'slm,ijkl,ijknm->ijksn', # tensor, 3x3 -> 3x3 - } - 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 # Nyquist freq=0 for even grid (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] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_si = np.arange(grid[0]//2+1)/geomdim[2] - - 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') - - 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, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0 - - curl_fourier = np.einsum(einsums[n],e,k_s,field_fourier)*TWOPIIMG - - return np.fft.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n]) - # -------------------------------------------------------------------- # MAIN @@ -60,8 +20,7 @@ def curlFFT(geomdim,field): parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """ Add column(s) containing curl of requested column(s). -Operates on periodic ordered three-dimensional data sets -of vector and tensor fields. +Operates on periodic ordered three-dimensional data sets of vector and tensor fields. """, version = scriptID) parser.add_option('-p','--pos','--periodiccellcenter', @@ -69,93 +28,30 @@ parser.add_option('-p','--pos','--periodiccellcenter', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') parser.add_option('-l','--label', - dest = 'data', + dest = 'labels', action = 'extend', metavar = '', help = 'label(s) of field values') parser.set_defaults(pos = 'pos', ) - (options,filenames) = parser.parse_args() - -if options.data is None: parser.error('no data column specified.') - -# --- define possible data types ------------------------------------------------------------------- - -datatypes = { - 3: {'name': 'vector', - 'shape': [3], - }, - 9: {'name': 'tensor', - 'shape': [3,3], - }, - } - -# --- loop over input files ------------------------------------------------------------------------ - if filenames == []: filenames = [None] +if options.labels is None: parser.error('no data column specified.') + for name in filenames: - try: table = damask.ASCIItable(name = name,buffered = False) - except: continue - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) -# --- interpret header ---------------------------------------------------------------------------- + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) + grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) - table.head_read() - - remarks = [] - errors = [] - active = [] - - coordDim = table.label_dimension(options.pos) - if coordDim != 3: - errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos)) - else: coordCol = table.label_index(options.pos) - - for me in options.data: - dim = table.label_dimension(me) - if dim in datatypes: - active.append(merge_dicts({'label':me},datatypes[dim])) - remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me)) - else: - remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \ - '"{}" not found...'.format(me) ) - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# ------------------------------------------ assemble header -------------------------------------- - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - for data in active: - table.labels_append(['{}_curlFFT({})'.format(i+1,data['label']) - for i in range(np.prod(np.array(data['shape'])))]) # extend ASCII header with new labels - table.head_write() - -# --------------- figure out size and grid --------------------------------------------------------- - - table.data_readArray() - grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) - -# ------------------------------------------ process value field ----------------------------------- - - stack = [table.data] - for data in active: - # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation - stack.append(curlFFT(size[::-1], - table.data[:,table.label_indexrange(data['label'])]. - reshape(grid[::-1].tolist()+data['shape']))) - -# ------------------------------------------ output result ----------------------------------------- - - if len(stack) > 1: table.data = np.hstack(tuple(stack)) - table.data_writeArray('%.12g') - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close input ASCII table (works for stdin) + for label in options.labels: + field = table.get_array(label) + shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor + field = table.get_array(label).reshape(np.append(grid[::-1],shape)) + table.add_array('curlFFT({})'.format(label), + damask.grid_filters.curl(size[::-1],field).reshape((-1,np.prod(shape))), + scriptID+' '+' '.join(sys.argv[1:])) + + table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 2d6af2036..562ab7532 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np @@ -12,52 +13,14 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -def merge_dicts(*dict_args): - """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts.""" - result = {} - for dictionary in dict_args: - result.update(dictionary) - return result - -def divFFT(geomdim,field): - """Calculate divergence of a vector or tensor field by transforming into Fourier space.""" - shapeFFT = np.array(np.shape(field))[0:3] - grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size - - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - - # differentiation in Fourier space - TWOPIIMG = 2.0j*np.pi - einsums = { - 3:'ijkl,ijkl->ijk', # vector, 3 -> 1 - 9:'ijkm,ijklm->ijkl', # tensor, 3x3 -> 3 - } - 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 # Nyquist freq=0 for even grid (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] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_si = np.arange(grid[0]//2+1)/geomdim[2] - - 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') - - div_fourier = np.einsum(einsums[n],k_s,field_fourier)*TWOPIIMG - - return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n//3]) - # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- -parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """ -Add column(s) containing curl of requested column(s). -Operates on periodic ordered three-dimensional data sets -of vector and tensor fields. +parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """ +Add column(s) containing divergence of requested column(s). +Operates on periodic ordered three-dimensional data sets of vector and tensor fields. """, version = scriptID) parser.add_option('-p','--pos','--periodiccellcenter', @@ -65,95 +28,30 @@ parser.add_option('-p','--pos','--periodiccellcenter', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') parser.add_option('-l','--label', - dest = 'data', + dest = 'labels', action = 'extend', metavar = '', help = 'label(s) of field values') parser.set_defaults(pos = 'pos', ) - (options,filenames) = parser.parse_args() - -if options.data is None: parser.error('no data column specified.') - -# --- define possible data types ------------------------------------------------------------------- - -datatypes = { - 3: {'name': 'vector', - 'shape': [3], - }, - 9: {'name': 'tensor', - 'shape': [3,3], - }, - } - -# --- loop over input files ------------------------------------------------------------------------ - if filenames == []: filenames = [None] +if options.labels is None: parser.error('no data column specified.') + for name in filenames: - try: table = damask.ASCIItable(name = name,buffered = False) - except: continue - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) -# --- interpret header ---------------------------------------------------------------------------- + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) + grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) - table.head_read() - - remarks = [] - errors = [] - active = [] - - coordDim = table.label_dimension(options.pos) - if coordDim != 3: - errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos)) - else: coordCol = table.label_index(options.pos) - - for me in options.data: - dim = table.label_dimension(me) - if dim in datatypes: - active.append(merge_dicts({'label':me},datatypes[dim])) - remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me)) - else: - remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \ - '"{}" not found...'.format(me) ) - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# ------------------------------------------ assemble header -------------------------------------- - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - for data in active: - table.labels_append(['divFFT({})'.format(data['label']) if data['shape'] == [3] \ - else '{}_divFFT({})'.format(i+1,data['label']) - for i in range(np.prod(np.array(data['shape']))//3)]) # extend ASCII header with new labels - table.head_write() - -# --------------- figure out size and grid --------------------------------------------------------- - - table.data_readArray() - - grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) - -# ------------------------------------------ process value field ----------------------------------- - - stack = [table.data] - for data in active: - # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation - stack.append(divFFT(size[::-1], - table.data[:,table.label_indexrange(data['label'])]. - reshape(grid[::-1].tolist()+data['shape']))) - -# ------------------------------------------ output result ----------------------------------------- - - if len(stack) > 1: table.data = np.hstack(tuple(stack)) - table.data_writeArray('%.12g') - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close input ASCII table (works for stdin) + for label in options.labels: + field = table.get_array(label) + shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor + field = table.get_array(label).reshape(np.append(grid[::-1],shape)) + table.add_array('divFFT({})'.format(label), + damask.grid_filters.divergence(size[::-1],field).reshape((-1,np.prod(shape)//3)), + scriptID+' '+' '.join(sys.argv[1:])) + + table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 252b72eb8..d6b537ddd 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np @@ -12,43 +13,6 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -def merge_dicts(*dict_args): - """Given any number of dicts, shallow copy and merge into a new dict, with precedence going to key value pairs in latter dicts.""" - result = {} - for dictionary in dict_args: - result.update(dictionary) - return result - -def gradFFT(geomdim,field): - """Calculate gradient of a vector or scalar field by transforming into Fourier space.""" - shapeFFT = np.array(np.shape(field))[0:3] - grid = np.array(np.shape(field)[2::-1]) - N = grid.prod() # field size - n = np.array(np.shape(field)[3:]).prod() # data size - - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=shapeFFT) - - # differentiation in Fourier space - TWOPIIMG = 2.0j*np.pi - einsums = { - 1:'ijkl,ijkm->ijkm', # scalar, 1 -> 3 - 3:'ijkl,ijkm->ijklm', # vector, 3 -> 3x3 - } - - 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 # Nyquist freq=0 for even grid (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] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_si = np.arange(grid[0]//2+1)/geomdim[2] - - 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') - grad_fourier = np.einsum(einsums[n],field_fourier,k_s)*TWOPIIMG - - return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n]) - # -------------------------------------------------------------------- # MAIN @@ -56,9 +20,7 @@ def gradFFT(geomdim,field): parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """ Add column(s) containing gradient of requested column(s). -Operates on periodic ordered three-dimensional data sets -of vector and scalar fields. - +Operates on periodic ordered three-dimensional data sets of scalar and vector fields. """, version = scriptID) parser.add_option('-p','--pos','--periodiccellcenter', @@ -66,7 +28,7 @@ parser.add_option('-p','--pos','--periodiccellcenter', type = 'string', metavar = 'string', help = 'label of coordinates [%default]') parser.add_option('-l','--label', - dest = 'data', + dest = 'labels', action = 'extend', metavar = '', help = 'label(s) of field values') @@ -74,85 +36,22 @@ parser.set_defaults(pos = 'pos', ) (options,filenames) = parser.parse_args() - -if options.data is None: parser.error('no data column specified.') - -# --- define possible data types ------------------------------------------------------------------- - -datatypes = { - 1: {'name': 'scalar', - 'shape': [1], - }, - 3: {'name': 'vector', - 'shape': [3], - }, - } - -# --- loop over input files ------------------------------------------------------------------------ - if filenames == []: filenames = [None] +if options.labels is None: parser.error('no data column specified.') + for name in filenames: - try: table = damask.ASCIItable(name = name,buffered = False) - except: continue - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) -# --- interpret header ---------------------------------------------------------------------------- + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) + grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) - table.head_read() - - remarks = [] - errors = [] - active = [] - - coordDim = table.label_dimension(options.pos) - if coordDim != 3: - errors.append('coordinates "{}" must be three-dimensional.'.format(options.pos)) - else: coordCol = table.label_index(options.pos) - - for me in options.data: - dim = table.label_dimension(me) - if dim in datatypes: - active.append(merge_dicts({'label':me},datatypes[dim])) - remarks.append('differentiating {} "{}"...'.format(datatypes[dim]['name'],me)) - else: - remarks.append('skipping "{}" of dimension {}...'.format(me,dim) if dim != -1 else \ - '"{}" not found...'.format(me) ) - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# ------------------------------------------ assemble header -------------------------------------- - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - for data in active: - table.labels_append(['{}_gradFFT({})'.format(i+1,data['label']) - for i in range(coordDim*np.prod(np.array(data['shape'])))]) # extend ASCII header with new labels - table.head_write() - -# --------------- figure out size and grid --------------------------------------------------------- - - table.data_readArray() - - grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)]) - -# ------------------------------------------ process value field ----------------------------------- - - stack = [table.data] - for data in active: - # we need to reverse order here, because x is fastest,ie rightmost, but leftmost in our x,y,z notation - stack.append(gradFFT(size[::-1], - table.data[:,table.label_indexrange(data['label'])]. - reshape(grid[::-1].tolist()+data['shape']))) - -# ------------------------------------------ output result ----------------------------------------- - - if len(stack) > 1: table.data = np.hstack(tuple(stack)) - table.data_writeArray('%.12g') - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close input ASCII table (works for stdin) + for label in options.labels: + field = table.get_array(label) + shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector + field = table.get_array(label).reshape(np.append(grid[::-1],shape)) + table.add_array('gradFFT({})'.format(label), + damask.grid_filters.gradient(size[::-1],field).reshape((-1,np.prod(shape)*3)), + scriptID+' '+' '.join(sys.argv[1:])) + + table.to_ASCII(sys.stdout if name is None else name) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 69ee85033..c7e96f468 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -1,12 +1,8 @@ import numpy as np - -def curl(size,field): - """Calculate curl of a vector or tensor field in Fourier space.""" +def __ks(size,field): + """Get differential operator.""" grid = np.array(np.shape(field)[0:3]) - n = np.array(np.shape(field)[3:]).prod() # data size - - field_fourier = np.fft.rfftn(field,axes=(0,1,2)) 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: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) @@ -17,64 +13,47 @@ def curl(size,field): k_si = np.arange(grid[2]//2+1)/size[2] 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') + return np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) + + +def curl(size,field): + """Calculate curl of a vector or tensor field in Fourier space.""" + n = np.prod(field.shape[3:]) + k_s = __ks(size,field) 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, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0 + field_fourier = np.fft.rfftn(field,axes=(0,1,2)) 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)) + return np.fft.irfftn(curl,axes=(0,1,2),s=field.shape[0:3]) def divergence(size,field): """Calculate divergence of a vector or tensor field in Fourier space.""" - grid = np.array(np.shape(field)[0:3]) - n = np.array(np.shape(field)[3:]).prod() # data size + n = np.prod(field.shape[3:]) + k_s = __ks(size,field) - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=grid) + 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 + np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*np.pi) # tensor, 3x3 -> 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: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_si = np.arange(grid[2]//2+1)/size[2] - - 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') - - 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(div_fourier,axes=(0,1,2),s=grid) + return np.fft.irfftn(divergence,axes=(0,1,2),s=field.shape[0:3]) def gradient(size,field): """Calculate gradient of a vector or scalar field in Fourier space.""" - grid = np.array(np.shape(field)[2::-1]) - n = np.array(np.shape(field)[3:]).prod() # data size - - field_fourier = np.fft.rfftn(field,axes=(0,1,2),s=grid) - - 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: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) - - k_si = np.arange(grid[2]//2+1)/size[2] - - 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') + n = np.prod(field.shape[3:]) + k_s = __ks(size,field) + 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 np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*np.pi) # vector, 3 -> 3x3 - return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=grid) + return np.fft.irfftn(gradient,axes=(0,1,2),s=field.shape[0:3]) #-------------------------------------------------------------------------------------------------- From f2e722ed2e92a9b61c9d18d8430619673ec1192c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Nov 2019 18:22:34 +0100 Subject: [PATCH 06/36] polishing --- python/damask/grid_filters.py | 64 ++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index c7e96f468..e49ff47a9 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -1,14 +1,14 @@ import numpy as np -def __ks(size,field): - """Get differential operator.""" +def __ks(size,field,first_order=False): + """Get wave numbers operator.""" grid = np.array(np.shape(field)[0: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: 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) k_sj = np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1]))/size[1] - if grid[1]%2 == 0: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) + if grid[1]%2 == 0 and first_order: k_sj[grid[1]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) k_si = np.arange(grid[2]//2+1)/size[2] @@ -19,7 +19,7 @@ def __ks(size,field): def curl(size,field): """Calculate curl of a vector or tensor field in Fourier space.""" n = np.prod(field.shape[3:]) - k_s = __ks(size,field) + k_s = __ks(size,field,True) e = np.zeros((3, 3, 3)) e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol @@ -35,7 +35,7 @@ def curl(size,field): def divergence(size,field): """Calculate divergence of a vector or tensor field in Fourier space.""" n = np.prod(field.shape[3:]) - k_s = __ks(size,field) + k_s = __ks(size,field,True) 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 @@ -47,7 +47,7 @@ def divergence(size,field): def gradient(size,field): """Calculate gradient of a vector or scalar field in Fourier space.""" n = np.prod(field.shape[3:]) - k_s = __ks(size,field) + k_s = __ks(size,field,True) 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 @@ -56,29 +56,37 @@ def gradient(size,field): return np.fft.irfftn(gradient,axes=(0,1,2),s=field.shape[0:3]) -#-------------------------------------------------------------------------------------------------- -def displacementFluctFFT(F,size): - """Calculate displacement field from deformation gradient field.""" - integrator = 0.5j * size / np.pi +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) - kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])), - np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])), - np.arange(grid[0]//2+1), - indexing = 'ij') - k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) - k_sSquared = np.einsum('...l,...l',k_s,k_s) - k_sSquared[0,0,0] = 1.0 # ignore global average frequency -#-------------------------------------------------------------------------------------------------- -# integration in Fourier space +def coord_cell(grid,size): + """Positions of cell centers (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]), + np.linspace(delta[0],size[0]-delta[0],grid[0]), + indexing = 'ij') - displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm', - np.fft.rfftn(F,axes=(0,1,2)), - k_s, - integrator, - ) / k_sSquared[...,np.newaxis] + return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) -#-------------------------------------------------------------------------------------------------- -# backtransformation to real space - return np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2)) +def displacement_fluct(size,F): + """Calculate displacement field from deformation gradient field.""" + integrator = 0.5j * size / np.pi + + k_s = __ks(size,F,False) + + displacement = -np.einsum('ijkml,ijkl,l->ijkm', + np.fft.rfftn(F,axes=(0,1,2)), + k_s, + integrator, + ) / k_sSquared[...,np.newaxis] + + return np.fft.irfftn(displacement,axes=(0,1,2)) From 62ca2952fce464016d721cf3ed2be5344d56a541 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 3 Dec 2019 11:27:14 +0100 Subject: [PATCH 07/36] polishing --- processing/post/addCurl.py | 2 +- processing/post/addDivergence.py | 2 +- processing/post/addGradient.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 2fcd107c0..b3dfedaf4 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -49,7 +49,7 @@ for name in filenames: for label in options.labels: field = table.get_array(label) shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor - field = table.get_array(label).reshape(np.append(grid[::-1],shape)) + field = field.reshape(np.append(grid[::-1],shape)) table.add_array('curlFFT({})'.format(label), damask.grid_filters.curl(size[::-1],field).reshape((-1,np.prod(shape))), scriptID+' '+' '.join(sys.argv[1:])) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 562ab7532..7ecaf10f0 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -49,7 +49,7 @@ for name in filenames: for label in options.labels: field = table.get_array(label) shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor - field = table.get_array(label).reshape(np.append(grid[::-1],shape)) + field = field.reshape(np.append(grid[::-1],shape)) table.add_array('divFFT({})'.format(label), damask.grid_filters.divergence(size[::-1],field).reshape((-1,np.prod(shape)//3)), scriptID+' '+' '.join(sys.argv[1:])) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index d6b537ddd..d3081db66 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -49,7 +49,7 @@ for name in filenames: for label in options.labels: field = table.get_array(label) shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector - field = table.get_array(label).reshape(np.append(grid[::-1],shape)) + field = field.reshape(np.append(grid[::-1],shape)) table.add_array('gradFFT({})'.format(label), damask.grid_filters.gradient(size[::-1],field).reshape((-1,np.prod(shape)*3)), scriptID+' '+' '.join(sys.argv[1:])) From e006e0ebec0a6e7a662f3f453040e2aa83d7f844 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 3 Dec 2019 18:59:59 +0100 Subject: [PATCH 08/36] functions for spatial coordinates on regular grids --- python/damask/grid_filters.py | 67 +++++++++++++++++++++++------------ 1 file changed, 45 insertions(+), 22 deletions(-) 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') From 9ad874339647b8271f5b66c5b9dc21554ee853fd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 3 Dec 2019 19:30:55 +0100 Subject: [PATCH 09/36] using central functionality --- processing/post/addDisplacement.py | 202 ++++------------------------- python/damask/grid_filters.py | 13 +- 2 files changed, 38 insertions(+), 177 deletions(-) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index 99d07fd18..ab25920b5 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -2,10 +2,10 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np -import scipy.ndimage import damask @@ -14,79 +14,6 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -#-------------------------------------------------------------------------------------------------- -def cell2node(cellData,grid): - - nodeData = 0.0 - datalen = np.array(cellData.shape[3:]).prod() - - for i in range(datalen): - 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 - ) # now averaged at cell origins - node = np.append(node,node[np.newaxis,0,:,:,...],axis=0) # wrap along z - node = np.append(node,node[:,0,np.newaxis,:,...],axis=1) # wrap along y - node = np.append(node,node[:,:,0,np.newaxis,...],axis=2) # wrap along x - - nodeData = node[...,np.newaxis] if i==0 else np.concatenate((nodeData,node[...,np.newaxis]),axis=-1) - - return nodeData - -#-------------------------------------------------------------------------------------------------- -def displacementAvgFFT(F,grid,size,nodal=False,transformed=False): - """Calculate average cell center (or nodal) displacement for deformation gradient field specified in each grid cell""" - if nodal: - 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') - else: - 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]), - np.linspace(delta[0],size[0]-delta[0],grid[0]), - indexing = 'ij') - - origCoords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) - - F_fourier = F if transformed else np.fft.rfftn(F,axes=(0,1,2)) # transform or use provided data - Favg = np.real(F_fourier[0,0,0,:,:])/grid.prod() # take zero freq for average - avgDisplacement = np.einsum('ml,ijkl->ijkm',Favg-np.eye(3),origCoords) # dX = Favg.X - - return avgDisplacement - -#-------------------------------------------------------------------------------------------------- -def displacementFluctFFT(F,grid,size,nodal=False,transformed=False): - """Calculate cell center (or nodal) displacement for deformation gradient field specified in each grid cell""" - integrator = 0.5j * size / np.pi - - kk, kj, ki = np.meshgrid(np.where(np.arange(grid[2])>grid[2]//2,np.arange(grid[2])-grid[2],np.arange(grid[2])), - np.where(np.arange(grid[1])>grid[1]//2,np.arange(grid[1])-grid[1],np.arange(grid[1])), - np.arange(grid[0]//2+1), - indexing = 'ij') - k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) - k_sSquared = np.einsum('...l,...l',k_s,k_s) - k_sSquared[0,0,0] = 1.0 # ignore global average frequency - -#-------------------------------------------------------------------------------------------------- -# integration in Fourier space - - displacement_fourier = -np.einsum('ijkml,ijkl,l->ijkm', - F if transformed else np.fft.rfftn(F,axes=(0,1,2)), - k_s, - integrator, - ) / k_sSquared[...,np.newaxis] - -#-------------------------------------------------------------------------------------------------- -# backtransformation to real space - - displacement = np.fft.irfftn(displacement_fourier,grid[::-1],axes=(0,1,2)) - - return cell2node(displacement,grid) if nodal else displacement - - # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -100,7 +27,7 @@ Outputs at cell centers or cell nodes (into separate file). parser.add_option('-f', '--defgrad', - dest = 'defgrad', + dest = 'f', metavar = 'string', help = 'label of deformation gradient [%default]') parser.add_option('-p', @@ -113,108 +40,35 @@ parser.add_option('--nodal', action = 'store_true', help = 'output nodal (instead of cell-centered) displacements') -parser.set_defaults(defgrad = 'f', - pos = 'pos', +parser.set_defaults(f = 'f', + pos = 'pos', ) (options,filenames) = parser.parse_args() -# --- loop over input files ------------------------------------------------------------------------- - -if filenames == []: filenames = [None] - for name in filenames: - outname = (os.path.splitext(name)[0] + - '_nodal' + - os.path.splitext(name)[1]) if (options.nodal and name) else None - try: table = damask.ASCIItable(name = name, - outname = outname, - buffered = False) - except: continue - damask.util.report(scriptName,'{}{}'.format(name if name else '', - ' --> {}'.format(outname) if outname else '')) + damask.util.report(scriptName,name) -# ------------------------------------------ read header ------------------------------------------ - - table.head_read() - -# ------------------------------------------ sanity checks ---------------------------------------- - - errors = [] - remarks = [] - - if table.label_dimension(options.defgrad) != 9: - errors.append('deformation gradient "{}" is not a 3x3 tensor.'.format(options.defgrad)) - - coordDim = table.label_dimension(options.pos) - if not 3 >= coordDim >= 1: - errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos)) - elif coordDim < 3: - remarks.append('appending {} dimension{} to coordinates "{}"...'.format(3-coordDim, - 's' if coordDim < 2 else '', - options.pos)) - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss=True) - continue - -# --------------- figure out size and grid --------------------------------------------------------- - - table.data_readArray([options.defgrad,options.pos]) - table.data_rewind() - - if len(table.data.shape) < 2: table.data.shape += (1,) # expand to 2D shape - if table.data[:,9:].shape[1] < 3: - table.data = np.hstack((table.data, - np.zeros((table.data.shape[0], - 3-table.data[:,9:].shape[1]),dtype='f'))) # fill coords up to 3D with zeros - - grid,size = damask.util.coordGridAndSize(table.data[:,9:12]) - N = grid.prod() - - if N != len(table.data): errors.append('data count {} does not match grid {}x{}x{}.'.format(N,*grid)) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# ------------------------------------------ process data ------------------------------------------ - - 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... - - fluctDisplacement = displacementFluctFFT(F_fourier,grid,size,options.nodal,transformed=True) - avgDisplacement = displacementAvgFFT (F_fourier,grid,size,options.nodal,transformed=True) - -# ------------------------------------------ assemble header --------------------------------------- - - if options.nodal: - table.info_clear() - table.labels_clear() - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - table.labels_append((['{}_pos' .format(i+1) for i in range(3)] if options.nodal else []) + - ['{}_avg({}).{}' .format(i+1,options.defgrad,options.pos) for i in range(3)] + - ['{}_fluct({}).{}'.format(i+1,options.defgrad,options.pos) for i in range(3)] ) - table.head_write() - -# ------------------------------------------ output data ------------------------------------------- - - Zrange = np.linspace(0,size[2],1+grid[2]) if options.nodal else range(grid[2]) - Yrange = np.linspace(0,size[1],1+grid[1]) if options.nodal else range(grid[1]) - Xrange = np.linspace(0,size[0],1+grid[0]) if options.nodal else range(grid[0]) - - 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(fluctDisplacement[i,j,k,:])) - table.data_write() - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close ASCII tables + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) + grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) + + F = table.get_array(options.f).reshape(np.append(grid[::-1],(3,3))) + if options.nodal: + table = damask.Table(damask.grid_filters.coord0_node(grid[::-1],size[::-1]).reshape((-1,3)), + {'pos':(3,)}) + table.add_array('avg({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_avg_node(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) + table.add_array('fluct({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_fluct_node(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) + table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt') + else: + table.add_array('avg({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_avg_cell(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) + table.add_array('fluct({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_fluct_cell(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) + + table.to_ASCII(sys.stdout if name is None else name) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index bf757bdb9..924242c58 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -102,14 +102,21 @@ def displacement_fluct_node(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)) + return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),coord0_node(F.shape[: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 + + 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') + +def node_2_cell(node_data): + """Interpolate nodal data to cell data.""" + 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,1)) + np.roll(node_data,1,(1,2)) + np.roll(node_data,1,(2,0)))*0.125 + + return c[:-1,:-1,:-1] From bc41bbbec52376d2a7b02449def6b0ffe6140c05 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 09:23:08 +0100 Subject: [PATCH 10/36] test coordinates-related functions --- python/tests/test_grid_filters.py | 42 +++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 python/tests/test_grid_filters.py diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py new file mode 100644 index 000000000..89b256dcc --- /dev/null +++ b/python/tests/test_grid_filters.py @@ -0,0 +1,42 @@ +import pytest +import numpy as np + +from damask import grid_filters + +class TestGridFilters: + + def test_coord0_cell(self): + size = np.random.random(3) + grid = np.random.randint(8,32,(3)) + coord = grid_filters.coord0_cell(grid,size) + assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid[::-1]) + (3,) + + def test_coord0_node(self): + size = np.random.random(3) + grid = np.random.randint(8,32,(3)) + coord = grid_filters.coord0_node(grid,size) + assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid[::-1]+1) + (3,) + + def test_coord0(self): + size = np.random.random(3) + grid = np.random.randint(8,32,(3)) + c = grid_filters.coord0_cell(grid+1,size+size/grid) + n = grid_filters.coord0_node(grid,size) + size/grid*.5 + assert np.allclose(c,n) + + @pytest.mark.parametrize('mode',[('cell'),('node')]) + def test_displacement_avg_vanishes(self,mode): + """Ensure that random fluctuations in F do not result in average displacement.""" + size = np.random.random(3) + grid = np.random.randint(8,32,(3)) + F = np.random.random(tuple(grid)+(3,3)) + F += np.eye(3) - np.average(F,axis=(0,1,2)) + assert np.allclose(eval('grid_filters.displacement_avg_{}(size,F)'.format(mode)),0.0) + + @pytest.mark.parametrize('mode',[('cell'),('node')]) + def test_displacement_fluct_vanishes(self,mode): + """Ensure that constant F does not result in fluctuating displacement.""" + size = np.random.random(3) + grid = np.random.randint(8,32,(3)) + F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) + assert np.allclose(eval('grid_filters.displacement_fluct_{}(size,F)'.format(mode)),0.0) From 381c95bd1e6c65600cb1d793942993c703085a9e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 4 Dec 2019 10:17:55 +0100 Subject: [PATCH 11/36] WIP: regrid functionality --- python/damask/grid_filters.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 924242c58..455e93aba 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -1,3 +1,4 @@ +from scipy import spatial import numpy as np def __ks(size,field,first_order=False): @@ -120,3 +121,18 @@ def node_2_cell(node_data): + 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] + +def regrid(size,F,new_grid): + """tbd.""" + c = coord0_cell(F.shape[:3][::-1],size) \ + + displacement_avg_cell(size,F) \ + + displacement_fluct_cell(size,F) + + outer = np.dot(np.average(F,axis=(0,1,2)),size) + for d in range(3): + c[np.where(c[:,:,:,d]<0)] += outer[d] + c[np.where(c[:,:,:,d]>outer[d])] -= outer[d] + + tree = spatial.cKDTree(c.reshape((-1,3)),boxsize=outer) + d,i = tree.query(coord0_cell(new_grid,outer)) + return i From e7358746dca1e8ab5f885374488121050b2ab469 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 5 Dec 2019 10:54:35 +0100 Subject: [PATCH 12/36] taking care of prospector complaints variables in 'eval' are hidden --- python/tests/test_grid_filters.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index 89b256dcc..e816e0380 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -27,7 +27,7 @@ class TestGridFilters: @pytest.mark.parametrize('mode',[('cell'),('node')]) def test_displacement_avg_vanishes(self,mode): """Ensure that random fluctuations in F do not result in average displacement.""" - size = np.random.random(3) + size = np.random.random(3) # noqa grid = np.random.randint(8,32,(3)) F = np.random.random(tuple(grid)+(3,3)) F += np.eye(3) - np.average(F,axis=(0,1,2)) @@ -36,7 +36,7 @@ class TestGridFilters: @pytest.mark.parametrize('mode',[('cell'),('node')]) def test_displacement_fluct_vanishes(self,mode): """Ensure that constant F does not result in fluctuating displacement.""" - size = np.random.random(3) + size = np.random.random(3) # noqa grid = np.random.randint(8,32,(3)) - F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) + F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) # noqa assert np.allclose(eval('grid_filters.displacement_fluct_{}(size,F)'.format(mode)),0.0) From f475d1a0d021912c080f652a72f9c32e50da8b95 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 5 Dec 2019 13:35:06 +0100 Subject: [PATCH 13/36] adjusted to changes in table class --- processing/post/addCurl.py | 8 ++++---- processing/post/addDisplacement.py | 27 +++++++++++++-------------- processing/post/addDivergence.py | 8 ++++---- processing/post/addGradient.py | 8 ++++---- 4 files changed, 25 insertions(+), 26 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index b3dfedaf4..6d4201cf3 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -47,11 +47,11 @@ for name in filenames: grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) for label in options.labels: - field = table.get_array(label) + field = table.get(label) shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor field = field.reshape(np.append(grid[::-1],shape)) - table.add_array('curlFFT({})'.format(label), - damask.grid_filters.curl(size[::-1],field).reshape((-1,np.prod(shape))), - scriptID+' '+' '.join(sys.argv[1:])) + table.add('curlFFT({})'.format(label), + damask.grid_filters.curl(size[::-1],field).reshape((-1,np.prod(shape))), + scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index ab25920b5..9c0e29fc5 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -52,23 +52,22 @@ for name in filenames: table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) - F = table.get_array(options.f).reshape(np.append(grid[::-1],(3,3))) + F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) if options.nodal: table = damask.Table(damask.grid_filters.coord0_node(grid[::-1],size[::-1]).reshape((-1,3)), {'pos':(3,)}) - table.add_array('avg({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_avg_node(size[::-1],F).reshape((-1,3)), - scriptID+' '+' '.join(sys.argv[1:])) - table.add_array('fluct({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_fluct_node(size[::-1],F).reshape((-1,3)), - scriptID+' '+' '.join(sys.argv[1:])) + table.add('avg({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_avg_node(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) + table.add('fluct({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_fluct_node(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt') else: - table.add_array('avg({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_avg_cell(size[::-1],F).reshape((-1,3)), - scriptID+' '+' '.join(sys.argv[1:])) - table.add_array('fluct({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_fluct_cell(size[::-1],F).reshape((-1,3)), - scriptID+' '+' '.join(sys.argv[1:])) - + table.add('avg({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_avg_cell(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) + table.add('fluct({}).{}'.format(options.f,options.pos), + damask.grid_filters.displacement_fluct_cell(size[::-1],F).reshape((-1,3)), + scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 7ecaf10f0..8c6f181f7 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -47,11 +47,11 @@ for name in filenames: grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) for label in options.labels: - field = table.get_array(label) + field = table.get(label) shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor field = field.reshape(np.append(grid[::-1],shape)) - table.add_array('divFFT({})'.format(label), - damask.grid_filters.divergence(size[::-1],field).reshape((-1,np.prod(shape)//3)), - scriptID+' '+' '.join(sys.argv[1:])) + table.add('divFFT({})'.format(label), + damask.grid_filters.divergence(size[::-1],field).reshape((-1,np.prod(shape)//3)), + scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index d3081db66..f7f800c24 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -47,11 +47,11 @@ for name in filenames: grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) for label in options.labels: - field = table.get_array(label) + field = table.get(label) shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector field = field.reshape(np.append(grid[::-1],shape)) - table.add_array('gradFFT({})'.format(label), - damask.grid_filters.gradient(size[::-1],field).reshape((-1,np.prod(shape)*3)), - scriptID+' '+' '.join(sys.argv[1:])) + table.add('gradFFT({})'.format(label), + damask.grid_filters.gradient(size[::-1],field).reshape((-1,np.prod(shape)*3)), + scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) From 4ddfd82304678f79d345d4535035a422254b2adf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 5 Dec 2019 18:32:21 +0100 Subject: [PATCH 14/36] better names --- python/damask/grid_filters.py | 26 +++++++++++++------------- python/tests/test_grid_filters.py | 16 ++++++++-------- 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 455e93aba..149d52194 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -57,7 +57,7 @@ def gradient(size,field): return np.fft.irfftn(gradient,axes=(0,1,2),s=field.shape[:3]) -def coord0_cell(grid,size): +def cell_coord0(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]), @@ -67,7 +67,7 @@ def coord0_cell(grid,size): return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) -def displacement_fluct_cell(size,F): +def cell_displacement_fluct(size,F): """Cell center displacement field from fluctuation part of the deformation gradient field.""" integrator = 0.5j*size/np.pi @@ -83,13 +83,13 @@ def displacement_fluct_cell(size,F): return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) -def displacement_avg_cell(size,F): +def cell_displacement_avg(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)) + return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),cell_coord0(F.shape[:3],size)) -def coord0_node(grid,size): +def node_coord0(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]), @@ -98,12 +98,12 @@ def coord0_node(grid,size): 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 node_displacement_fluct(size,F): + return cell_2_node(cell_displacement_fluct(size,F)) -def displacement_avg_node(size,F): +def node_displacement_avg(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[:3],size)) + return np.einsum('ml,ijkl->ijkm',F_avg-np.eye(3),node_coord0(F.shape[:3],size)) def cell_2_node(cell_data): @@ -124,9 +124,9 @@ def node_2_cell(node_data): def regrid(size,F,new_grid): """tbd.""" - c = coord0_cell(F.shape[:3][::-1],size) \ - + displacement_avg_cell(size,F) \ - + displacement_fluct_cell(size,F) + c = cell_coord0(F.shape[:3][::-1],size) \ + + cell_displacement_avg(size,F) \ + + cell_displacement_fluct(size,F) outer = np.dot(np.average(F,axis=(0,1,2)),size) for d in range(3): @@ -134,5 +134,5 @@ def regrid(size,F,new_grid): c[np.where(c[:,:,:,d]>outer[d])] -= outer[d] tree = spatial.cKDTree(c.reshape((-1,3)),boxsize=outer) - d,i = tree.query(coord0_cell(new_grid,outer)) + d,i = tree.query(cell_coord0(new_grid,outer)) return i diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index e816e0380..4dff41542 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -5,23 +5,23 @@ from damask import grid_filters class TestGridFilters: - def test_coord0_cell(self): + def test_cell_coord0(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) - coord = grid_filters.coord0_cell(grid,size) + coord = grid_filters.cell_coord0(grid,size) assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid[::-1]) + (3,) - def test_coord0_node(self): + def test_node_coord0(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) - coord = grid_filters.coord0_node(grid,size) + coord = grid_filters.node_coord0(grid,size) assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid[::-1]+1) + (3,) def test_coord0(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) - c = grid_filters.coord0_cell(grid+1,size+size/grid) - n = grid_filters.coord0_node(grid,size) + size/grid*.5 + c = grid_filters.cell_coord0(grid+1,size+size/grid) + n = grid_filters.node_coord0(grid,size) + size/grid*.5 assert np.allclose(c,n) @pytest.mark.parametrize('mode',[('cell'),('node')]) @@ -31,7 +31,7 @@ class TestGridFilters: grid = np.random.randint(8,32,(3)) F = np.random.random(tuple(grid)+(3,3)) F += np.eye(3) - np.average(F,axis=(0,1,2)) - assert np.allclose(eval('grid_filters.displacement_avg_{}(size,F)'.format(mode)),0.0) + assert np.allclose(eval('grid_filters.{}_displacement_avg(size,F)'.format(mode)),0.0) @pytest.mark.parametrize('mode',[('cell'),('node')]) def test_displacement_fluct_vanishes(self,mode): @@ -39,4 +39,4 @@ class TestGridFilters: size = np.random.random(3) # noqa grid = np.random.randint(8,32,(3)) F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) # noqa - assert np.allclose(eval('grid_filters.displacement_fluct_{}(size,F)'.format(mode)),0.0) + assert np.allclose(eval('grid_filters.{}_displacement_fluct(size,F)'.format(mode)),0.0) From f2ac87eb2f60e27881217f1586508c7df6904055 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 6 Dec 2019 04:22:18 +0100 Subject: [PATCH 15/36] follow changes in Table class --- processing/post/addCurl.py | 2 +- processing/post/addDisplacement.py | 2 +- processing/post/addDivergence.py | 2 +- processing/post/addGradient.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 6d4201cf3..010d24935 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) + grid,size = damask.util.coordGridAndSize(table.get(options.pos)) for label in options.labels: field = table.get(label) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index 9c0e29fc5..f85c2a914 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -50,7 +50,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) + grid,size = damask.util.coordGridAndSize(table.get(options.pos)) F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) if options.nodal: diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 8c6f181f7..d2d68ede6 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) + grid,size = damask.util.coordGridAndSize(table.get(options.pos)) for label in options.labels: field = table.get(label) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index f7f800c24..ed9397f02 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get_array(options.pos)) + grid,size = damask.util.coordGridAndSize(table.get(options.pos)) for label in options.labels: field = table.get(label) From 3aab154cdbef90ae4ba43657b336290a240b9d03 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 7 Dec 2019 10:55:06 +0100 Subject: [PATCH 16/36] include test for addDisplacement --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 524e86c11..806e1b32a 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 524e86c117d816e3bd873eed7663e258a6f2e139 +Subproject commit 806e1b32a17bf26c987f417116476a295d3936cd From e283acd6069d7f729f38bc835f2c33a92fd7986e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 7 Dec 2019 20:07:46 +0100 Subject: [PATCH 17/36] correct way of importing for newer python versions --- processing/post/addCalculation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/processing/post/addCalculation.py b/processing/post/addCalculation.py index db0428753..6effc2489 100755 --- a/processing/post/addCalculation.py +++ b/processing/post/addCalculation.py @@ -4,7 +4,7 @@ import os import sys from optparse import OptionParser import re -import collections +from collections.abc import Iterable import math # noqa import scipy # noqa @@ -18,7 +18,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) def listify(x): - return x if isinstance(x, collections.Iterable) else [x] + return x if isinstance(x, Iterable) else [x] # -------------------------------------------------------------------- From c6c77b64d2b7d5347e90c5d334fec321cb7bee9f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 7 Dec 2019 20:08:31 +0100 Subject: [PATCH 18/36] following renames in grid_filter --- processing/post/addDisplacement.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index f85c2a914..4f0323504 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -54,20 +54,20 @@ for name in filenames: F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) if options.nodal: - table = damask.Table(damask.grid_filters.coord0_node(grid[::-1],size[::-1]).reshape((-1,3)), + table = damask.Table(damask.grid_filters.node_coord0(grid[::-1],size[::-1]).reshape((-1,3)), {'pos':(3,)}) table.add('avg({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_avg_node(size[::-1],F).reshape((-1,3)), + damask.grid_filters.node_displacement_avg(size[::-1],F).reshape((-1,3)), scriptID+' '+' '.join(sys.argv[1:])) table.add('fluct({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_fluct_node(size[::-1],F).reshape((-1,3)), + damask.grid_filters.node_displacement_fluct(size[::-1],F).reshape((-1,3)), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt') else: table.add('avg({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_avg_cell(size[::-1],F).reshape((-1,3)), + damask.grid_filters.cell_displacement_avg(size[::-1],F).reshape((-1,3)), scriptID+' '+' '.join(sys.argv[1:])) table.add('fluct({}).{}'.format(options.f,options.pos), - damask.grid_filters.displacement_fluct_cell(size[::-1],F).reshape((-1,3)), + damask.grid_filters.cell_displacement_fluct(size[::-1],F).reshape((-1,3)), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) From b92cfbbd5bce653ec0bf0a5a0518f1c87bbfb60d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 7 Dec 2019 20:15:50 +0100 Subject: [PATCH 19/36] do not use bare except --- processing/post/addCalculation.py | 7 ++++--- processing/post/addCumulative.py | 6 +++--- processing/post/addOrientations.py | 7 ++++--- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/processing/post/addCalculation.py b/processing/post/addCalculation.py index 6effc2489..b1eed3c6d 100755 --- a/processing/post/addCalculation.py +++ b/processing/post/addCalculation.py @@ -65,9 +65,10 @@ for i in range(len(options.formulas)): if filenames == []: filenames = [None] for name in filenames: - try: table = damask.ASCIItable(name = name, - buffered = False) - except: continue + try: + table = damask.ASCIItable(name = name, buffered = False) + except IOError: + continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------- diff --git a/processing/post/addCumulative.py b/processing/post/addCumulative.py index b81a9d14f..c94737b94 100755 --- a/processing/post/addCumulative.py +++ b/processing/post/addCumulative.py @@ -41,9 +41,9 @@ if filenames == []: filenames = [None] for name in filenames: try: - table = damask.ASCIItable(name = name, - buffered = False) - except IOError: continue + table = damask.ASCIItable(name = name, buffered = False) + except IOError: + continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 31ce6aeb3..2c46ee5ee 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -125,9 +125,10 @@ R = damask.Rotation.fromAxisAngle(np.array(options.labrotation),options.degrees, if filenames == []: filenames = [None] for name in filenames: - try: table = damask.ASCIItable(name = name, - buffered = False) - except Exception: continue + try: + table = damask.ASCIItable(name = name, buffered = False) + except IOError: + continue damask.util.report(scriptName,name) # ------------------------------------------ read header ------------------------------------------ From ba69f5a631735075b372618503b0c31eed069a8e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 7 Dec 2019 22:33:31 +0100 Subject: [PATCH 20/36] polishing --- python/damask/grid_filters.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 149d52194..006167ccc 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -99,9 +99,11 @@ def node_coord0(grid,size): return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) def node_displacement_fluct(size,F): + """Nodal displacement field from fluctuation part of the deformation gradient field.""" return cell_2_node(cell_displacement_fluct(size,F)) def node_displacement_avg(size,F): + """Nodal 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),node_coord0(F.shape[:3],size)) @@ -134,5 +136,4 @@ def regrid(size,F,new_grid): c[np.where(c[:,:,:,d]>outer[d])] -= outer[d] tree = spatial.cKDTree(c.reshape((-1,3)),boxsize=outer) - d,i = tree.query(cell_coord0(new_grid,outer)) - return i + return tree.query(cell_coord0(new_grid,outer))[1] From 9dc726ff53bf43e8e33f95d263e2e33bdbfe0524 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 09:17:57 +0100 Subject: [PATCH 21/36] polishing --- python/damask/geom.py | 3 +++ python/damask/table.py | 1 - 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/python/damask/geom.py b/python/damask/geom.py index 32ea2ed89..63ce20115 100644 --- a/python/damask/geom.py +++ b/python/damask/geom.py @@ -205,6 +205,9 @@ class Geom(): else: self.homogenization = homogenization + @property + def grid(self): + return self.get_grid() def get_microstructure(self): """Return the microstructure representation.""" diff --git a/python/damask/table.py b/python/damask/table.py index 6181fdb1f..d063599ab 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -97,7 +97,6 @@ class Table(): @property def labels(self): - """Return the labels of all columns.""" return list(self.shapes.keys()) From 12564557e65e83e6510ed626119d9a023220f0ee Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 09:18:15 +0100 Subject: [PATCH 22/36] using central functionality --- processing/pre/seeds_fromGeom.py | 106 +++++++++---------------------- 1 file changed, 29 insertions(+), 77 deletions(-) diff --git a/processing/pre/seeds_fromGeom.py b/processing/pre/seeds_fromGeom.py index 889ef6146..c543cb878 100755 --- a/processing/pre/seeds_fromGeom.py +++ b/processing/pre/seeds_fromGeom.py @@ -1,9 +1,12 @@ #!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- -import os,sys -import numpy as np +import os +import sys +from io import StringIO from optparse import OptionParser + +import numpy as np + import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] @@ -29,88 +32,37 @@ parser.add_option('-b', action = 'extend', metavar = '', dest = 'blacklist', help = 'blacklist of grain IDs') -parser.add_option('-p', - '--pos', '--seedposition', - dest = 'pos', - type = 'string', metavar = 'string', - help = 'label of coordinates [%default]') parser.set_defaults(whitelist = [], blacklist = [], - pos = 'pos', ) (options,filenames) = parser.parse_args() - -options.whitelist = list(map(int,options.whitelist)) -options.blacklist = list(map(int,options.blacklist)) - -# --- loop over output files ------------------------------------------------------------------------- - if filenames == []: filenames = [None] +options.whitelist = [int(i) for i in options.whitelist] +options.blacklist = [int(i) for i in options.blacklist] + for name in filenames: - try: table = damask.ASCIItable(name = name, - outname = os.path.splitext(name)[0]+'.seeds' if name else name, - buffered = False, - labeled = False) - except: continue - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) + + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) + microstructure = geom.get_microstructure().reshape((-1,1),order='F') -# --- interpret header ---------------------------------------------------------------------------- + mask = np.logical_and(np.in1d(microstructure,options.whitelist,invert=False) if options.whitelist else \ + np.full(geom.grid.prod(),True,dtype=bool), + np.in1d(microstructure,options.blacklist,invert=True) if options.blacklist else \ + np.full(geom.grid.prod(),True,dtype=bool)) + + seeds = np.concatenate((damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape((-1,3)), + microstructure), + axis=1)[mask] + + comments = [scriptID + ' ' + ' '.join(sys.argv[1:]), + "grid\ta {}\tb {}\tc {}".format(*geom.grid), + "size\tx {}\ty {}\tz {}".format(*geom.size), + "origin\tx {}\ty {}\tz {}".format(*geom.origin), + "homogenization\t{}".format(geom.homogenization)] - table.head_read() - info,extra_header = table.head_getGeom() - damask.util.report_geom(info) - - errors = [] - if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') - if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.') - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# --- read data ------------------------------------------------------------------------------------ - - microstructure = table.microstructure_read(info['grid']) # read (linear) microstructure - -# --- generate grid -------------------------------------------------------------------------------- - - x = (0.5 + np.arange(info['grid'][0],dtype=float))/info['grid'][0]*info['size'][0]+info['origin'][0] - y = (0.5 + np.arange(info['grid'][1],dtype=float))/info['grid'][1]*info['size'][1]+info['origin'][1] - z = (0.5 + np.arange(info['grid'][2],dtype=float))/info['grid'][2]*info['size'][2]+info['origin'][2] - - xx = np.tile( x, info['grid'][1]* info['grid'][2]) - yy = np.tile(np.repeat(y,info['grid'][0] ),info['grid'][2]) - zz = np.repeat(z,info['grid'][0]*info['grid'][1]) - - mask = np.logical_and(np.in1d(microstructure,options.whitelist,invert=False) if options.whitelist != [] - else np.full_like(microstructure,True,dtype=bool), - np.in1d(microstructure,options.blacklist,invert=True ) if options.blacklist != [] - else np.full_like(microstructure,True,dtype=bool)) - -# ------------------------------------------ assemble header --------------------------------------- - - table.info_clear() - table.info_append(extra_header+[ - scriptID + ' ' + ' '.join(sys.argv[1:]), - "grid\ta {}\tb {}\tc {}".format(*info['grid']), - "size\tx {}\ty {}\tz {}".format(*info['size']), - "origin\tx {}\ty {}\tz {}".format(*info['origin']), - "homogenization\t{}".format(info['homogenization']), - "microstructures\t{}".format(info['microstructures']), - ]) - table.labels_clear() - table.labels_append(['{dim}_{label}'.format(dim = 1+i,label = options.pos) for i in range(3)]+['microstructure']) - table.head_write() - table.output_flush() - -# --- write seeds information ------------------------------------------------------------ - - table.data = np.squeeze(np.dstack((xx,yy,zz,microstructure)))[mask] - table.data_writeArray() - -# ------------------------------------------ finalize output --------------------------------------- - - table.close() + table = damask.Table(seeds,{'pos':(3,),'microstructure':(1,)},comments) + table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds') From 871ff4c218363649be2ad307e2144698c6cb5d5b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 09:31:56 +0100 Subject: [PATCH 23/36] use geom class --- processing/pre/mentat_spectralBox.py | 101 ++++++++++----------------- 1 file changed, 35 insertions(+), 66 deletions(-) diff --git a/processing/pre/mentat_spectralBox.py b/processing/pre/mentat_spectralBox.py index 89f4a7a43..f62622a6b 100755 --- a/processing/pre/mentat_spectralBox.py +++ b/processing/pre/mentat_spectralBox.py @@ -1,9 +1,11 @@ #!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- -import os,sys -import numpy as np +import os +import sys from optparse import OptionParser + +import numpy as np + import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] @@ -191,78 +193,45 @@ parser.add_option('-p', '--port', dest = 'port', type = 'int', metavar = 'int', help = 'Mentat connection port [%default]') -parser.add_option('--homogenization', - dest = 'homogenization', - type = 'int', metavar = 'int', - help = 'homogenization index to be used [auto]') -parser.set_defaults(port = None, - homogenization = None, -) +parser.set_defaults(port = None, + ) (options, filenames) = parser.parse_args() -if options.port: - try: - import py_mentat - except: - parser.error('no valid Mentat release found.') +if options.port is not None: + try: + import py_mentat + except ImportError: + parser.error('no valid Mentat release found.') # --- loop over input files ------------------------------------------------------------------------ if filenames == []: filenames = [None] for name in filenames: - try: - table = damask.ASCIItable(name = name, - outname = os.path.splitext(name)[0]+'.proc' if name else name, - buffered = False, labeled = False) - except: continue - damask.util.report(scriptName,name) - -# --- interpret header ---------------------------------------------------------------------------- - - table.head_read() - info,extra_header = table.head_getGeom() - if options.homogenization: info['homogenization'] = options.homogenization + damask.util.report(scriptName,name) - damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))), - 'size x y z: %s'%(' x '.join(map(str,info['size']))), - 'origin x y z: %s'%(' : '.join(map(str,info['origin']))), - 'homogenization: %i'%info['homogenization'], - 'microstructures: %i'%info['microstructures'], - ]) + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) + microstructure = geom.get_microstructure().flatten(order='F') - errors = [] - if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') - if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.') - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# --- read data ------------------------------------------------------------------------------------ - - microstructure = table.microstructure_read(info['grid']).reshape(info['grid'].prod(),order='F') # read microstructure - - cmds = [\ - init(), - mesh(info['grid'],info['size']), - material(), - geometry(), - initial_conditions(info['homogenization'],microstructure), - '*identify_sets', - '*show_model', - '*redraw', - '*draw_automatic', - ] - - outputLocals = {} - if options.port: - py_mentat.py_connect('',options.port) - output(cmds,outputLocals,'Mentat') - py_mentat.py_disconnect() - else: - output(cmds,outputLocals,table.__IO__['out']) # bad hack into internals of table class... - - table.close() + cmds = [\ + init(), + mesh(geom.grid,geom.size), + material(), + geometry(), + initial_conditions(geom.homogenization,microstructure), + '*identify_sets', + '*show_model', + '*redraw', + '*draw_automatic', + ] + + outputLocals = {} + if options.port: + py_mentat.py_connect('',options.port) + output(cmds,outputLocals,'Mentat') + py_mentat.py_disconnect() + else: + with sys.stdout if name is None else open(os.path.splitext(name)[0]+'.proc','w') as f: + output(cmds,outputLocals,f) From 0292e8fcc78ff7c87036885192ed38a7e61758bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 11:00:38 +0100 Subject: [PATCH 24/36] preparing transition to Geom and Table classes --- processing/pre/seeds_fromPokes.py | 100 ++++++++++-------------------- 1 file changed, 33 insertions(+), 67 deletions(-) diff --git a/processing/pre/seeds_fromPokes.py b/processing/pre/seeds_fromPokes.py index 08e600ffe..4a6edc94d 100755 --- a/processing/pre/seeds_fromPokes.py +++ b/processing/pre/seeds_fromPokes.py @@ -1,11 +1,13 @@ #!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- -import os,math,sys -import numpy as np -import damask +import os +import sys from optparse import OptionParser +import numpy as np + +import damask + scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) @@ -35,109 +37,73 @@ parser.add_option('-y', action = 'store_true', dest = 'y', help = 'poke 45 deg along y') -parser.add_option('-p','--position', - dest = 'position', - type = 'string', metavar = 'string', - help = 'column label for coordinates [%default]') parser.set_defaults(x = False, y = False, box = [0.0,1.0,0.0,1.0,0.0,1.0], N = 16, - position = 'pos', ) (options,filenames) = parser.parse_args() +if filenames == []: filenames = [None] options.box = np.array(options.box).reshape(3,2) -# --- loop over output files ------------------------------------------------------------------------- - -if filenames == []: filenames = [None] - for name in filenames: - try: - table = damask.ASCIItable(name = name, - outname = os.path.splitext(name)[-2]+'_poked_{}.seeds'.format(options.N) if name else name, - buffered = False, labeled = False) - except: continue + table = damask.ASCIItable(name = name, + outname = os.path.splitext(name)[-2]+'_poked_{}.seeds'.format(options.N) if name else name, + buffered = False, labeled = False) damask.util.report(scriptName,name) -# --- interpret header ---------------------------------------------------------------------------- table.head_read() info,extra_header = table.head_getGeom() - - damask.util.croak(['grid a b c: %s'%(' x '.join(map(str,info['grid']))), - 'size x y z: %s'%(' x '.join(map(str,info['size']))), - 'origin x y z: %s'%(' : '.join(map(str,info['origin']))), - 'homogenization: %i'%info['homogenization'], - 'microstructures: %i'%info['microstructures'], - ]) - - errors = [] - if np.any(info['grid'] < 1): errors.append('invalid grid a b c.') - if np.any(info['size'] <= 0.0): errors.append('invalid size x y z.') - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# --- read data ------------------------------------------------------------------------------------ - - microstructure = table.microstructure_read(info['grid']).reshape(info['grid'],order='F') # read microstructure - -# --- do work ------------------------------------------------------------------------------------ - newInfo = { - 'microstructures': 0, - } - offset = (np.amin(options.box, axis=1)*info['grid']/info['size']).astype(int) - box = np.amax(options.box, axis=1) - np.amin(options.box, axis=1) + grid = info['grid'] + size = info['size'] - Nx = int(options.N/math.sqrt(options.N*info['size'][1]*box[1]/info['size'][0]/box[0])) - Ny = int(options.N/math.sqrt(options.N*info['size'][0]*box[0]/info['size'][1]/box[1])) - Nz = int(box[2]*info['grid'][2]) + microstructure = table.microstructure_read(grid).reshape(grid) # read microstructure + + offset =(np.amin(options.box, axis=1)*grid/size).astype(int) + box = np.amax(options.box, axis=1) \ + - np.amin(options.box, axis=1) + + Nx = int(options.N/np.sqrt(options.N*size[1]*box[1]/size[0]/box[0])) + Ny = int(options.N/np.sqrt(options.N*size[0]*box[0]/size[1]/box[1])) + Nz = int(box[2]*grid[2]) damask.util.croak('poking {} x {} x {} in box {} {} {}...'.format(Nx,Ny,Nz,*box)) seeds = np.zeros((Nx*Ny*Nz,4),'d') - grid = np.zeros(3,'i') + g = np.zeros(3,'i') n = 0 for i in range(Nx): for j in range(Ny): - grid[0] = round((i+0.5)*box[0]*info['grid'][0]/Nx-0.5)+offset[0] - grid[1] = round((j+0.5)*box[1]*info['grid'][1]/Ny-0.5)+offset[1] + g[0] = round((i+0.5)*box[0]*grid[0]/Nx-0.5)+offset[0] + g[1] = round((j+0.5)*box[1]*grid[1]/Ny-0.5)+offset[1] for k in range(Nz): - grid[2] = k + offset[2] - grid %= info['grid'] - seeds[n,0:3] = (0.5+grid)/info['grid'] # normalize coordinates to box - seeds[n, 3] = microstructure[grid[0],grid[1],grid[2]] - if options.x: grid[0] += 1 - if options.y: grid[1] += 1 + g[2] = k + offset[2] + g %= grid + seeds[n,0:3] = (g+0.5)/grid # normalize coordinates to box + seeds[n, 3] = microstructure[g[2],g[1],g[0]] + if options.x: g[0] += 1 + if options.y: g[1] += 1 n += 1 - newInfo['microstructures'] = len(np.unique(seeds[:,3])) - -# --- report --------------------------------------------------------------------------------------- - if (newInfo['microstructures'] != info['microstructures']): - damask.util.croak('--> microstructures: %i'%newInfo['microstructures']) - # ------------------------------------------ assemble header --------------------------------------- table.info_clear() table.info_append(extra_header+[ scriptID + ' ' + ' '.join(sys.argv[1:]), "poking\ta {}\tb {}\tc {}".format(Nx,Ny,Nz), - "grid\ta {}\tb {}\tc {}".format(*info['grid']), - "size\tx {}\ty {}\tz {}".format(*info['size']), + "grid\ta {}\tb {}\tc {}".format(*grid), + "size\tx {}\ty {}\tz {}".format(*size), "origin\tx {}\ty {}\tz {}".format(*info['origin']), "homogenization\t{}".format(info['homogenization']), - "microstructures\t{}".format(newInfo['microstructures']), ]) table.labels_clear() - table.labels_append(['{dim}_{label}'.format(dim = 1+i,label = options.position) for i in range(3)]+['microstructure']) + table.labels_append(['{dim}_{label}'.format(dim = 1+i,label = 'pos') for i in range(3)]+['microstructure']) table.head_write() table.output_flush() From f19694f734c401e4076af1294951d2ef140c45d2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 11:20:47 +0100 Subject: [PATCH 25/36] starting to use central functionality --- processing/pre/seeds_fromPokes.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/processing/pre/seeds_fromPokes.py b/processing/pre/seeds_fromPokes.py index 4a6edc94d..2c359b0aa 100755 --- a/processing/pre/seeds_fromPokes.py +++ b/processing/pre/seeds_fromPokes.py @@ -50,11 +50,12 @@ if filenames == []: filenames = [None] options.box = np.array(options.box).reshape(3,2) for name in filenames: + damask.util.report(scriptName,name) + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) + table = damask.ASCIItable(name = name, outname = os.path.splitext(name)[-2]+'_poked_{}.seeds'.format(options.N) if name else name, buffered = False, labeled = False) - damask.util.report(scriptName,name) - table.head_read() info,extra_header = table.head_getGeom() @@ -62,15 +63,13 @@ for name in filenames: grid = info['grid'] size = info['size'] - microstructure = table.microstructure_read(grid).reshape(grid) # read microstructure - - offset =(np.amin(options.box, axis=1)*grid/size).astype(int) + offset =(np.amin(options.box, axis=1)*geom.grid/geom.size).astype(int) box = np.amax(options.box, axis=1) \ - np.amin(options.box, axis=1) - Nx = int(options.N/np.sqrt(options.N*size[1]*box[1]/size[0]/box[0])) - Ny = int(options.N/np.sqrt(options.N*size[0]*box[0]/size[1]/box[1])) - Nz = int(box[2]*grid[2]) + Nx = int(options.N/np.sqrt(options.N*geom.size[1]*box[1]/geom.size[0]/box[0])) + Ny = int(options.N/np.sqrt(options.N*geom.size[0]*box[0]/geom.size[1]/box[1])) + Nz = int(box[2]*geom.grid[2]) damask.util.croak('poking {} x {} x {} in box {} {} {}...'.format(Nx,Ny,Nz,*box)) @@ -86,7 +85,7 @@ for name in filenames: g[2] = k + offset[2] g %= grid seeds[n,0:3] = (g+0.5)/grid # normalize coordinates to box - seeds[n, 3] = microstructure[g[2],g[1],g[0]] + seeds[n, 3] = geom.microstructure[g[0],g[1],g[2]] if options.x: g[0] += 1 if options.y: g[1] += 1 n += 1 @@ -94,12 +93,12 @@ for name in filenames: # ------------------------------------------ assemble header --------------------------------------- table.info_clear() - table.info_append(extra_header+[ + table.info_append(geom.comments+[ scriptID + ' ' + ' '.join(sys.argv[1:]), "poking\ta {}\tb {}\tc {}".format(Nx,Ny,Nz), - "grid\ta {}\tb {}\tc {}".format(*grid), - "size\tx {}\ty {}\tz {}".format(*size), - "origin\tx {}\ty {}\tz {}".format(*info['origin']), + "grid\ta {}\tb {}\tc {}".format(*geom.grid), + "size\tx {}\ty {}\tz {}".format(*geom.size), + "origin\tx {}\ty {}\tz {}".format(*geom.origin), "homogenization\t{}".format(info['homogenization']), ]) table.labels_clear() From 75e93d9f0c96d88fb2b204fcd14f000cd7954cf2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 11:25:33 +0100 Subject: [PATCH 26/36] relying on central functionality improves readability --- processing/pre/seeds_fromGeom.py | 6 +- processing/pre/seeds_fromPokes.py | 104 ++++++++++++------------------ 2 files changed, 45 insertions(+), 65 deletions(-) diff --git a/processing/pre/seeds_fromGeom.py b/processing/pre/seeds_fromGeom.py index c543cb878..2118f049d 100755 --- a/processing/pre/seeds_fromGeom.py +++ b/processing/pre/seeds_fromGeom.py @@ -58,11 +58,13 @@ for name in filenames: microstructure), axis=1)[mask] - comments = [scriptID + ' ' + ' '.join(sys.argv[1:]), + comments = geom.comments \ + + [scriptID + ' ' + ' '.join(sys.argv[1:]), "grid\ta {}\tb {}\tc {}".format(*geom.grid), "size\tx {}\ty {}\tz {}".format(*geom.size), "origin\tx {}\ty {}\tz {}".format(*geom.origin), "homogenization\t{}".format(geom.homogenization)] table = damask.Table(seeds,{'pos':(3,),'microstructure':(1,)},comments) - table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds') + table.to_ASCII(sys.stdout if name is None else \ + os.path.splitext(name)[0]+'.seeds') diff --git a/processing/pre/seeds_fromPokes.py b/processing/pre/seeds_fromPokes.py index 2c359b0aa..ef7da63bf 100755 --- a/processing/pre/seeds_fromPokes.py +++ b/processing/pre/seeds_fromPokes.py @@ -50,67 +50,45 @@ if filenames == []: filenames = [None] options.box = np.array(options.box).reshape(3,2) for name in filenames: - damask.util.report(scriptName,name) - geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) - - table = damask.ASCIItable(name = name, - outname = os.path.splitext(name)[-2]+'_poked_{}.seeds'.format(options.N) if name else name, - buffered = False, labeled = False) - - table.head_read() - info,extra_header = table.head_getGeom() - - grid = info['grid'] - size = info['size'] - - offset =(np.amin(options.box, axis=1)*geom.grid/geom.size).astype(int) - box = np.amax(options.box, axis=1) \ - - np.amin(options.box, axis=1) - - Nx = int(options.N/np.sqrt(options.N*geom.size[1]*box[1]/geom.size[0]/box[0])) - Ny = int(options.N/np.sqrt(options.N*geom.size[0]*box[0]/geom.size[1]/box[1])) - Nz = int(box[2]*geom.grid[2]) - - damask.util.croak('poking {} x {} x {} in box {} {} {}...'.format(Nx,Ny,Nz,*box)) - - seeds = np.zeros((Nx*Ny*Nz,4),'d') - g = np.zeros(3,'i') - - n = 0 - for i in range(Nx): - for j in range(Ny): - g[0] = round((i+0.5)*box[0]*grid[0]/Nx-0.5)+offset[0] - g[1] = round((j+0.5)*box[1]*grid[1]/Ny-0.5)+offset[1] - for k in range(Nz): - g[2] = k + offset[2] - g %= grid - seeds[n,0:3] = (g+0.5)/grid # normalize coordinates to box - seeds[n, 3] = geom.microstructure[g[0],g[1],g[2]] - if options.x: g[0] += 1 - if options.y: g[1] += 1 - n += 1 - - -# ------------------------------------------ assemble header --------------------------------------- - table.info_clear() - table.info_append(geom.comments+[ - scriptID + ' ' + ' '.join(sys.argv[1:]), - "poking\ta {}\tb {}\tc {}".format(Nx,Ny,Nz), - "grid\ta {}\tb {}\tc {}".format(*geom.grid), - "size\tx {}\ty {}\tz {}".format(*geom.size), - "origin\tx {}\ty {}\tz {}".format(*geom.origin), - "homogenization\t{}".format(info['homogenization']), - ]) - table.labels_clear() - table.labels_append(['{dim}_{label}'.format(dim = 1+i,label = 'pos') for i in range(3)]+['microstructure']) - table.head_write() - table.output_flush() - -# --- write seeds information ------------------------------------------------------------ - - table.data = seeds - table.data_writeArray() + damask.util.report(scriptName,name) + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) -# --- output finalization -------------------------------------------------------------------------- - - table.close() # close ASCII table + offset =(np.amin(options.box, axis=1)*geom.grid/geom.size).astype(int) + box = np.amax(options.box, axis=1) \ + - np.amin(options.box, axis=1) + + Nx = int(options.N/np.sqrt(options.N*geom.size[1]*box[1]/geom.size[0]/box[0])) + Ny = int(options.N/np.sqrt(options.N*geom.size[0]*box[0]/geom.size[1]/box[1])) + Nz = int(box[2]*geom.grid[2]) + + damask.util.croak('poking {} x {} x {} in box {} {} {}...'.format(Nx,Ny,Nz,*box)) + + seeds = np.zeros((Nx*Ny*Nz,4),'d') + g = np.zeros(3,'i') + + n = 0 + for i in range(Nx): + for j in range(Ny): + g[0] = round((i+0.5)*box[0]*geom.grid[0]/Nx-0.5)+offset[0] + g[1] = round((j+0.5)*box[1]*geom.grid[1]/Ny-0.5)+offset[1] + for k in range(Nz): + g[2] = k + offset[2] + g %= geom.grid + seeds[n,0:3] = (g+0.5)/geom.grid # normalize coordinates to box + seeds[n, 3] = geom.microstructure[g[0],g[1],g[2]] + if options.x: g[0] += 1 + if options.y: g[1] += 1 + n += 1 + + + comments = geom.comments \ + + [scriptID + ' ' + ' '.join(sys.argv[1:]), + "poking\ta {}\tb {}\tc {}".format(Nx,Ny,Nz), + "grid\ta {}\tb {}\tc {}".format(*geom.grid), + "size\tx {}\ty {}\tz {}".format(*geom.size), + "origin\tx {}\ty {}\tz {}".format(*geom.origin), + "homogenization\t{}".format(geom.homogenization)] + + table = damask.Table(seeds,{'pos':(3,),'microstructure':(1,)},comments) + table.to_ASCII(sys.stdout if name is None else \ + os.path.splitext(name)[0]+'_poked_{}.seeds'.format(options.N)) From 8d0c4310cf54acf3182258095ad68527de845498 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 17:57:02 +0100 Subject: [PATCH 27/36] improvements to grid generation - handling of origin - inverse functions: calculate grid,size,origin from regular coordinates (cell or node). should replace corresponding functionality in util --- python/damask/grid_filters.py | 66 ++++++++++++++++++++++++++++++----- 1 file changed, 57 insertions(+), 9 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 006167ccc..cb593f449 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -57,12 +57,13 @@ def gradient(size,field): return np.fft.irfftn(gradient,axes=(0,1,2),s=field.shape[:3]) -def cell_coord0(grid,size): +def cell_coord0(grid,size,origin=np.zeros(3)): """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]), - np.linspace(delta[0],size[0]-delta[0],grid[0]), + start = origin + size/grid*.5 + end = origin - size/grid*.5 + size + x, y, z = np.meshgrid(np.linspace(start[2],end[2],grid[2]), + np.linspace(start[1],end[1],grid[1]), + np.linspace(start[0],end[0],grid[0]), indexing = 'ij') return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) @@ -88,12 +89,37 @@ def cell_displacement_avg(size,F): 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],size)) +def cell_coord0_2_DNA(coord0,ordered=False): + coords = [np.unique(coord0[:,i]) for i in range(3)] + mincorner = np.array(list(map(min,coords))) + maxcorner = np.array(list(map(max,coords))) + grid = np.array(list(map(len,coords)),'i') + size = grid/np.maximum(grid-1,1) * (maxcorner-mincorner) + delta = size/grid + origin = mincorner - delta*.5 + + if grid.prod() != len(coord0): + raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) -def node_coord0(grid,size): + start = origin + delta*.5 + end = origin + size -delta*.5 + + if not np.allclose(coords[0],np.linspace(start[0],end[0],grid[0])) and \ + np.allclose(coords[1],np.linspace(start[1],end[1],grid[1])) and \ + np.allclose(coords[2],np.linspace(start[2],end[2],grid[2])): + raise ValueError('Regular grid spacing violated.') + + if ordered: + if not np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): + raise ValueError('Input data is not a regular grid.') + return (grid,size,origin) + + +def node_coord0(grid,size,origin=np.zeros(3)): """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]), + x, y, z = np.meshgrid(np.linspace(origin[2],size[2]+origin[2],1+grid[2]), + np.linspace(origin[1],size[1]+origin[1],1+grid[1]), + np.linspace(origin[0],size[0]+origin[0],1+grid[0]), indexing = 'ij') return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) @@ -124,6 +150,28 @@ def node_2_cell(node_data): return c[:-1,:-1,:-1] +def node_coord0_2_DNA(coord0,ordered=False): + coords = [np.unique(coord0[:,i]) for i in range(3)] + mincorner = np.array(list(map(min,coords))) + maxcorner = np.array(list(map(max,coords))) + grid = np.array(list(map(len,coords)),'i') - 1 + size = maxcorner-mincorner + origin = mincorner + + if (grid+1).prod() != len(coord0): + raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) + + if not np.allclose(coords[0],np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ + np.allclose(coords[1],np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ + np.allclose(coords[2],np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): + raise ValueError('Regular grid spacing violated.') + + if ordered: + if not np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): + raise ValueError('Input data is not a regular grid.') + return (grid,size,origin) + + def regrid(size,F,new_grid): """tbd.""" c = cell_coord0(F.shape[:3][::-1],size) \ From 828e82605ed401be736af10cd033b8ef474cd653 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 18:13:20 +0100 Subject: [PATCH 28/36] ensure that data is correctly ordered --- processing/post/addCurl.py | 2 +- processing/post/addDisplacement.py | 2 +- processing/post/addDivergence.py | 2 +- processing/post/addGradient.py | 2 +- python/damask/grid_filters.py | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 010d24935..5db47ce04 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get(options.pos)) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) for label in options.labels: field = table.get(label) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index 4f0323504..735e6d875 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -50,7 +50,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get(options.pos)) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) if options.nodal: diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index d2d68ede6..ea4b134a4 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get(options.pos)) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) for label in options.labels: field = table.get(label) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index ed9397f02..75109a3e0 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size = damask.util.coordGridAndSize(table.get(options.pos)) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) for label in options.labels: field = table.get(label) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index cb593f449..96ae226dd 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -102,7 +102,7 @@ def cell_coord0_2_DNA(coord0,ordered=False): raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) start = origin + delta*.5 - end = origin + size -delta*.5 + end = origin - delta*.5 + size if not np.allclose(coords[0],np.linspace(start[0],end[0],grid[0])) and \ np.allclose(coords[1],np.linspace(start[1],end[1],grid[1])) and \ From b0689340d060567d21b8b2b4ca0cdb3fc7bd0e60 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 18:20:28 +0100 Subject: [PATCH 29/36] updated tests --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 806e1b32a..b580fd4e9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 806e1b32a17bf26c987f417116476a295d3936cd +Subproject commit b580fd4e992c2ed457f16d65bcba2e6d099dc29d From 3d09a82f41898d5deb96a8a35967a7f253726f2a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 18:21:16 +0100 Subject: [PATCH 30/36] fixing prospector complaints --- processing/pre/mentat_spectralBox.py | 3 +-- processing/pre/seeds_fromPokes.py | 1 + 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/processing/pre/mentat_spectralBox.py b/processing/pre/mentat_spectralBox.py index f62622a6b..a61bef57a 100755 --- a/processing/pre/mentat_spectralBox.py +++ b/processing/pre/mentat_spectralBox.py @@ -2,10 +2,9 @@ import os import sys +from io import StringIO from optparse import OptionParser -import numpy as np - import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] diff --git a/processing/pre/seeds_fromPokes.py b/processing/pre/seeds_fromPokes.py index ef7da63bf..1436841d0 100755 --- a/processing/pre/seeds_fromPokes.py +++ b/processing/pre/seeds_fromPokes.py @@ -2,6 +2,7 @@ import os import sys +from io import StringIO from optparse import OptionParser import numpy as np From 7dc128ad123afb3b2d9e4b608e7855d5d724903b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 18:33:43 +0100 Subject: [PATCH 31/36] polishing --- processing/post/addCurl.py | 2 +- processing/post/addDisplacement.py | 2 +- processing/post/addDivergence.py | 2 +- processing/post/addGradient.py | 2 +- python/damask/grid_filters.py | 36 ++++++++++++++++++++++++------ 5 files changed, 33 insertions(+), 11 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 5db47ce04..25639dc7c 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos)) for label in options.labels: field = table.get(label) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index 735e6d875..59630a6c6 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -50,7 +50,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos)) F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) if options.nodal: diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index ea4b134a4..585ebb5a5 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos)) for label in options.labels: field = table.get(label) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 75109a3e0..54b80ed26 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -44,7 +44,7 @@ for name in filenames: damask.util.report(scriptName,name) table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) - grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos),True) + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.get(options.pos)) for label in options.labels: field = table.get(label) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 96ae226dd..b064d9a2d 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -89,7 +89,18 @@ def cell_displacement_avg(size,F): 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],size)) -def cell_coord0_2_DNA(coord0,ordered=False): +def cell_coord0_2_DNA(coord0,ordered=True): + """ + Return grid 'DNA', i.e. grid, size, and origin from array of cell positions. + + Parameters + ---------- + coord0 : numpy.ndarray + array of undeformed cell coordinates + ordered : bool, optional + expect coord0 data to be ordered (x fast, z slow). + + """ coords = [np.unique(coord0[:,i]) for i in range(3)] mincorner = np.array(list(map(min,coords))) maxcorner = np.array(list(map(max,coords))) @@ -109,9 +120,9 @@ def cell_coord0_2_DNA(coord0,ordered=False): np.allclose(coords[2],np.linspace(start[2],end[2],grid[2])): raise ValueError('Regular grid spacing violated.') - if ordered: - if not np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): - raise ValueError('Input data is not a regular grid.') + if ordered and not np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): + raise ValueError('Input data is not a regular grid.') + return (grid,size,origin) @@ -151,6 +162,17 @@ def node_2_cell(node_data): return c[:-1,:-1,:-1] def node_coord0_2_DNA(coord0,ordered=False): + """ + Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions. + + Parameters + ---------- + coord0 : numpy.ndarray + array of undeformed nodal coordinates + ordered : bool, optional + expect coord0 data to be ordered (x fast, z slow). + + """ coords = [np.unique(coord0[:,i]) for i in range(3)] mincorner = np.array(list(map(min,coords))) maxcorner = np.array(list(map(max,coords))) @@ -166,9 +188,9 @@ def node_coord0_2_DNA(coord0,ordered=False): np.allclose(coords[2],np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): raise ValueError('Regular grid spacing violated.') - if ordered: - if not np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): - raise ValueError('Input data is not a regular grid.') + if ordered and not np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): + raise ValueError('Input data is not a regular grid.') + return (grid,size,origin) From 4e2e7d02f6c10ad80a2ce9f4fa49257e723ceec8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 18:54:41 +0100 Subject: [PATCH 32/36] more sensible interface --- python/damask/grid_filters.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index b064d9a2d..9e1170ec4 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -1,9 +1,8 @@ from scipy import spatial import numpy as np -def __ks(size,field,first_order=False): +def __ks(size,grid,first_order=False): """Get wave numbers operator.""" - 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) @@ -20,7 +19,7 @@ def __ks(size,field,first_order=False): def curl(size,field): """Calculate curl of a vector or tensor field in Fourier space.""" n = np.prod(field.shape[3:]) - k_s = __ks(size,field,True) + k_s = __ks(size,field.shape[:3],True) e = np.zeros((3, 3, 3)) e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol @@ -36,7 +35,7 @@ def curl(size,field): def divergence(size,field): """Calculate divergence of a vector or tensor field in Fourier space.""" n = np.prod(field.shape[3:]) - k_s = __ks(size,field,True) + k_s = __ks(size,field.shape[:3],True) 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 @@ -48,7 +47,7 @@ def divergence(size,field): def gradient(size,field): """Calculate gradient of a vector or scalar field in Fourier space.""" n = np.prod(field.shape[3:]) - k_s = __ks(size,field,True) + k_s = __ks(size,field.shape[:3],True) 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 @@ -72,7 +71,7 @@ def cell_displacement_fluct(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 = __ks(size,F.shape[:3],False) k_s_squared = np.einsum('...l,...l',k_s,k_s) k_s_squared[0,0,0] = 1.0 From 21431295fbf8231e3d014f87a13deb7ac1b46a0d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 20:20:13 +0100 Subject: [PATCH 33/36] documenting --- python/damask/grid_filters.py | 113 ++++++++++++++++++++++++++++++---- 1 file changed, 102 insertions(+), 11 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 9e1170ec4..a1e1ff06d 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -2,8 +2,15 @@ from scipy import spatial import numpy as np def __ks(size,grid,first_order=False): - """Get wave numbers operator.""" + """ + Get wave numbers operator. + Parameters + ---------- + size : numpy.ndarray + 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] if grid[0]%2 == 0 and first_order: k_sk[grid[0]//2] = 0 # Nyquist freq=0 for even grid (Johnson, MIT, 2011) @@ -17,7 +24,15 @@ def __ks(size,grid,first_order=False): def curl(size,field): - """Calculate curl of a vector or tensor field in Fourier space.""" + """ + Calculate curl of a vector or tensor field in Fourier space. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + + """ n = np.prod(field.shape[3:]) k_s = __ks(size,field.shape[:3],True) @@ -33,7 +48,15 @@ def curl(size,field): def divergence(size,field): - """Calculate divergence of a vector or tensor field in Fourier space.""" + """ + Calculate divergence of a vector or tensor field in Fourier space. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + + """ n = np.prod(field.shape[3:]) k_s = __ks(size,field.shape[:3],True) @@ -45,7 +68,15 @@ def divergence(size,field): def gradient(size,field): - """Calculate gradient of a vector or scalar field in Fourier space.""" + """ + Calculate gradient of a vector or scalar field in Fourier space. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + + """ n = np.prod(field.shape[3:]) k_s = __ks(size,field.shape[:3],True) @@ -57,7 +88,17 @@ def gradient(size,field): def cell_coord0(grid,size,origin=np.zeros(3)): - """Cell center positions (undeformed).""" + """ + Cell center positions (undeformed). + + Parameters + ---------- + grid : numpy.ndarray + number of grid points. + size : numpy.ndarray + physical size of the periodic field. + + """ start = origin + size/grid*.5 end = origin - size/grid*.5 + size x, y, z = np.meshgrid(np.linspace(start[2],end[2],grid[2]), @@ -68,7 +109,17 @@ def cell_coord0(grid,size,origin=np.zeros(3)): return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) def cell_displacement_fluct(size,F): - """Cell center displacement field from fluctuation part of the deformation gradient field.""" + """ + Cell center displacement field from fluctuation part of the deformation gradient field. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + F : numpy.ndarray + deformation gradient field. + + """ integrator = 0.5j*size/np.pi k_s = __ks(size,F.shape[:3],False) @@ -84,7 +135,17 @@ def cell_displacement_fluct(size,F): return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) 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. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + F : numpy.ndarray + deformation gradient field. + + """ 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],size)) @@ -95,7 +156,7 @@ def cell_coord0_2_DNA(coord0,ordered=True): Parameters ---------- coord0 : numpy.ndarray - array of undeformed cell coordinates + array of undeformed cell coordinates. ordered : bool, optional expect coord0 data to be ordered (x fast, z slow). @@ -126,7 +187,17 @@ def cell_coord0_2_DNA(coord0,ordered=True): def node_coord0(grid,size,origin=np.zeros(3)): - """Nodal positions (undeformed).""" + """ + Nodal positions (undeformed). + + Parameters + ---------- + grid : numpy.ndarray + number of grid points. + size : numpy.ndarray + physical size of the periodic field. + + """ x, y, z = np.meshgrid(np.linspace(origin[2],size[2]+origin[2],1+grid[2]), np.linspace(origin[1],size[1]+origin[1],1+grid[1]), np.linspace(origin[0],size[0]+origin[0],1+grid[0]), @@ -135,11 +206,31 @@ def node_coord0(grid,size,origin=np.zeros(3)): return np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) def node_displacement_fluct(size,F): - """Nodal displacement field from fluctuation part of the deformation gradient field.""" + """ + Nodal displacement field from fluctuation part of the deformation gradient field. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + F : numpy.ndarray + deformation gradient field. + + """ return cell_2_node(cell_displacement_fluct(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. + + Parameters + ---------- + size : numpy.ndarray + physical size of the periodic field. + F : numpy.ndarray + deformation gradient field. + + """ 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],size)) From b56864552f92ce306b23fcbf26a839bdf3c547d8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 21:05:34 +0100 Subject: [PATCH 34/36] testing forward <-> backward conversion --- python/tests/test_grid_filters.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index 4dff41542..b23fad549 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -24,6 +24,18 @@ class TestGridFilters: n = grid_filters.node_coord0(grid,size) + size/grid*.5 assert np.allclose(c,n) + @pytest.mark.parametrize('mode',[('cell'),('node')]) + def test_grid_DNA(self,mode): + """Ensure that xx_coord0_2_DNA is the inverse of xx_coord0.""" + grid = np.random.randint(8,32,(3)) + size = np.random.random(3) + origin = np.random.random(3) + + coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) # noqa + _grid,_size,_origin = eval('grid_filters.{}_coord0_2_DNA(coord0.reshape((-1,3)))'.format(mode)) + assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin) + + @pytest.mark.parametrize('mode',[('cell'),('node')]) def test_displacement_avg_vanishes(self,mode): """Ensure that random fluctuations in F do not result in average displacement.""" From 0bf22fd03c9e5512d61b4954df9eac0733a0ecd6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 22:35:39 +0100 Subject: [PATCH 35/36] using central functionality --- processing/pre/geom_toTable.py | 45 ++++++++++++---------------------- 1 file changed, 15 insertions(+), 30 deletions(-) diff --git a/processing/pre/geom_toTable.py b/processing/pre/geom_toTable.py index ed33d9c85..3d955ad53 100755 --- a/processing/pre/geom_toTable.py +++ b/processing/pre/geom_toTable.py @@ -2,10 +2,8 @@ import os import sys -from optparse import OptionParser from io import StringIO - -import numpy as np +from optparse import OptionParser import damask @@ -24,38 +22,25 @@ Translate geom description into ASCIItable containing position and microstructur """, version = scriptID) (options, filenames) = parser.parse_args() - - if filenames == []: filenames = [None] for name in filenames: - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) - geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) + damask.util.croak(geom) - damask.util.croak(geom) + coord0 = damask.grid_filters.cell_coord0(geom.grid,geom.size,geom.origin).reshape((-1,3),order='F') -# --- generate grid -------------------------------------------------------------------------------- + comments = geom.comments \ + + [scriptID + ' ' + ' '.join(sys.argv[1:]), + "grid\ta {}\tb {}\tc {}".format(*geom.grid), + "size\tx {}\ty {}\tz {}".format(*geom.size), + "origin\tx {}\ty {}\tz {}".format(*geom.origin), + "homogenization\t{}".format(geom.homogenization)] - grid = geom.get_grid() - size = geom.get_size() - origin = geom.get_origin() + table = damask.Table(coord0,{'pos':(3,)},comments) + table.add('microstructure',geom.microstructure.reshape((-1,1))) - x = (0.5 + np.arange(grid[0],dtype=float))/grid[0]*size[0]+origin[0] - y = (0.5 + np.arange(grid[1],dtype=float))/grid[1]*size[1]+origin[1] - z = (0.5 + np.arange(grid[2],dtype=float))/grid[2]*size[2]+origin[2] - - xx = np.tile( x, grid[1]* grid[2]) - yy = np.tile(np.repeat(y,grid[0] ),grid[2]) - zz = np.repeat(z,grid[0]*grid[1]) - -# --- create ASCII table -------------------------------------------------------------------------- - - table = damask.ASCIItable(outname = os.path.splitext(name)[0]+'.txt' if name else name) - table.info_append(geom.get_comments() + [scriptID + '\t' + ' '.join(sys.argv[1:])]) - table.labels_append(['{}_{}'.format(1+i,'pos') for i in range(3)]+['microstructure']) - table.head_write() - table.output_flush() - table.data = np.squeeze(np.dstack((xx,yy,zz,geom.microstructure.flatten('F'))),axis=0) - table.data_writeArray() - table.close() + table.to_ASCII(sys.stdout if name is None else \ + os.path.splitext(name)[0]+'.txt') From 1955d09d3f92a3f6ddcd3956689567942eb1ef8c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 8 Dec 2019 22:35:58 +0100 Subject: [PATCH 36/36] modernizing --- processing/pre/geom_fromTable.py | 24 +----------------------- 1 file changed, 1 insertion(+), 23 deletions(-) diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index f513c4834..40fd17437 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -78,36 +78,15 @@ for name in filenames: table = damask.ASCIItable(name = name,readonly=True) table.head_read() # read ASCII header info -# ------------------------------------------ sanity checks --------------------------------------- - - coordDim = table.label_dimension(options.pos) - - errors = [] - if not 3 >= coordDim >= 2: - errors.append('coordinates "{}" need to have two or three dimensions.'.format(options.pos)) - if not np.all(table.label_dimension(label) == dim): - errors.append('input "{}" needs to have dimension {}.'.format(label,dim)) - if options.phase and table.label_dimension(options.phase) != 1: - errors.append('phase column "{}" is not scalar.'.format(options.phase)) - - if errors != []: - damask.util.croak(errors) - continue table.data_readArray([options.pos] \ + (label if isinstance(label, list) else [label]) \ + ([options.phase] if options.phase else [])) - if coordDim == 2: - table.data = np.insert(table.data,2,np.zeros(len(table.data)),axis=1) # add zero z coordinate for two-dimensional input if options.phase is None: table.data = np.column_stack((table.data,np.ones(len(table.data)))) # add single phase if no phase column given - grid,size = damask.util.coordGridAndSize(table.data[:,0:3]) - coords = [np.unique(table.data[:,i]) for i in range(3)] - mincorner = np.array(list(map(min,coords))) - origin = mincorner - 0.5*size/grid # shift from cell center to corner - + grid,size,origin = damask.grid_filters.cell_coord0_2_DNA(table.data[:,0:3]) indices = np.lexsort((table.data[:,0],table.data[:,1],table.data[:,2])) # indices of position when sorting x fast, z slow microstructure = np.empty(grid,dtype = int) # initialize empty microstructure @@ -142,7 +121,6 @@ for name in filenames: config_header += [''] for i,data in enumerate(unique): config_header += ['[Grain{}]'.format(i+1), - 'crystallite 1', '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(int(data[4]),i+1), ]