diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 87c557721..9a01e8e62 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -6,7 +6,7 @@ name = 'damask' with open(_os.path.join(_os.path.dirname(__file__),'VERSION')) as _f: version = _re.sub(r'^v','',_f.readline().strip()) -# classes +# make classes directly accessible as damask.Class from ._environment import Environment # noqa from ._table import Table # noqa from ._vtk import VTK # noqa diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 7288be6cb..e2c2428fe 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -291,7 +291,7 @@ class Geom: comments = [] for i,line in enumerate(content[:header_length]): - items = line.lower().strip().split() + items = line.split('#')[0].lower().strip().split() key = items[0] if items else '' if key == 'grid': grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']]) @@ -307,7 +307,7 @@ class Geom: microstructure = np.empty(grid.prod()) # initialize as flat array i = 0 for line in content[header_length:]: - items = line.split() + items = line.split('#')[0].split() if len(items) == 3: if items[1].lower() == 'of': items = np.ones(int(items[0]))*float(items[2]) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 35876fc99..8106ed905 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -1,5 +1,5 @@ -from scipy import spatial -import numpy as np +from scipy import spatial as _spatial +import numpy as _np def _ks(size,grid,first_order=False): """ @@ -11,16 +11,16 @@ def _ks(size,grid,first_order=False): physical size of the periodic field. """ - k_sk = np.where(np.arange(grid[0])>grid[0]//2,np.arange(grid[0])-grid[0],np.arange(grid[0]))/size[0] + k_sk = _np.where(_np.arange(grid[0])>grid[0]//2,_np.arange(grid[0])-grid[0],_np.arange(grid[0]))/size[0] 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] + 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 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] + k_si = _np.arange(grid[2]//2+1)/size[2] - kk, kj, ki = np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') - return np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) + kk, kj, ki = _np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') + return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) def curl(size,field): @@ -33,18 +33,18 @@ def curl(size,field): physical size of the periodic field. """ - n = np.prod(field.shape[3:]) + n = _np.prod(field.shape[3:]) k_s = _ks(size,field.shape[:3],True) - e = np.zeros((3, 3, 3)) + e = _np.zeros((3, 3, 3)) e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol e[0, 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 + 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),s=field.shape[:3]) + return _np.fft.irfftn(curl_,axes=(0,1,2),s=field.shape[:3]) def divergence(size,field): @@ -57,14 +57,14 @@ def divergence(size,field): physical size of the periodic field. """ - n = np.prod(field.shape[3:]) + n = _np.prod(field.shape[3:]) k_s = _ks(size,field.shape[:3],True) - field_fourier = np.fft.rfftn(field,axes=(0,1,2)) - div_ = (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 + field_fourier = _np.fft.rfftn(field,axes=(0,1,2)) + div_ = (_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_,axes=(0,1,2),s=field.shape[:3]) + return _np.fft.irfftn(div_,axes=(0,1,2),s=field.shape[:3]) def gradient(size,field): @@ -77,17 +77,17 @@ def gradient(size,field): physical size of the periodic field. """ - n = np.prod(field.shape[3:]) + n = _np.prod(field.shape[3:]) k_s = _ks(size,field.shape[:3],True) - field_fourier = np.fft.rfftn(field,axes=(0,1,2)) - grad_ = (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 + field_fourier = _np.fft.rfftn(field,axes=(0,1,2)) + grad_ = (_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_,axes=(0,1,2),s=field.shape[:3]) + return _np.fft.irfftn(grad_,axes=(0,1,2),s=field.shape[:3]) -def cell_coord0(grid,size,origin=np.zeros(3)): +def cell_coord0(grid,size,origin=_np.zeros(3)): """ Cell center positions (undeformed). @@ -103,7 +103,7 @@ def cell_coord0(grid,size,origin=np.zeros(3)): """ start = origin + size/grid*.5 end = origin + size - size/grid*.5 - return np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T + return _np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T def cell_displacement_fluct(size,F): @@ -118,19 +118,19 @@ def cell_displacement_fluct(size,F): deformation gradient field. """ - integrator = 0.5j*size/np.pi + integrator = 0.5j*size/_np.pi k_s = _ks(size,F.shape[:3],False) - k_s_squared = np.einsum('...l,...l',k_s,k_s) + k_s_squared = _np.einsum('...l,...l',k_s,k_s) k_s_squared[0,0,0] = 1.0 - displacement = -np.einsum('ijkml,ijkl,l->ijkm', - np.fft.rfftn(F,axes=(0,1,2)), + displacement = -_np.einsum('ijkml,ijkl,l->ijkm', + _np.fft.rfftn(F,axes=(0,1,2)), k_s, integrator, - ) / k_s_squared[...,np.newaxis] + ) / k_s_squared[...,_np.newaxis] - return np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) + return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) def cell_displacement_avg(size,F): @@ -145,8 +145,8 @@ def cell_displacement_avg(size,F): 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][::-1],size)) + F_avg = _np.average(F,axis=(0,1,2)) + return _np.einsum('ml,ijkl->ijkm',F_avg-_np.eye(3),cell_coord0(F.shape[:3][::-1],size)) def cell_displacement(size,F): @@ -164,7 +164,7 @@ def cell_displacement(size,F): return cell_displacement_avg(size,F) + cell_displacement_fluct(size,F) -def cell_coord(size,F,origin=np.zeros(3)): +def cell_coord(size,F,origin=_np.zeros(3)): """ Cell center positions. @@ -193,17 +193,17 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): 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))) - grid = np.array(list(map(len,coords)),'i') - size = grid/np.maximum(grid-1,1) * (maxcorner-mincorner) + 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 # 1D/2D: size/origin combination undefined, set origin to 0.0 - size [np.where(grid==1)] = origin[np.where(grid==1)]*2. - origin[np.where(grid==1)] = 0.0 + size [_np.where(grid==1)] = origin[_np.where(grid==1)]*2. + origin[_np.where(grid==1)] = 0.0 if grid.prod() != len(coord0): raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) @@ -211,13 +211,13 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): start = origin + 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 \ - np.allclose(coords[2],np.linspace(start[2],end[2],grid[2])): + 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 and 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('I_nput data is not a regular grid.') return (grid,size,origin) @@ -235,7 +235,7 @@ def coord0_check(coord0): cell_coord0_gridSizeOrigin(coord0,ordered=True) -def node_coord0(grid,size,origin=np.zeros(3)): +def node_coord0(grid,size,origin=_np.zeros(3)): """ Nodal positions (undeformed). @@ -249,7 +249,7 @@ def node_coord0(grid,size,origin=np.zeros(3)): physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. """ - return np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j, + return _np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j, origin[1]:size[1]+origin[1]:(grid[1]+1)*1j, origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T @@ -281,8 +281,8 @@ def node_displacement_avg(size,F): 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][::-1],size)) + F_avg = _np.average(F,axis=(0,1,2)) + return _np.einsum('ml,ijkl->ijkm',F_avg-_np.eye(3),node_coord0(F.shape[:3][::-1],size)) def node_displacement(size,F): @@ -300,7 +300,7 @@ def node_displacement(size,F): return node_displacement_avg(size,F) + node_displacement_fluct(size,F) -def node_coord(size,F,origin=np.zeros(3)): +def node_coord(size,F,origin=_np.zeros(3)): """ Nodal positions. @@ -319,18 +319,18 @@ def node_coord(size,F,origin=np.zeros(3)): def cell_2_node(cell_data): """Interpolate periodic 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 + 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') + 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 periodic 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 + 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] @@ -347,23 +347,23 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False): 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))) - grid = np.array(list(map(len,coords)),'i') - 1 + 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)): + 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 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.') + if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): + raise ValueError('I_nput data is not a regular grid.') return (grid,size,origin) @@ -386,10 +386,10 @@ def regrid(size,F,new_grid): + cell_displacement_avg(size,F) \ + cell_displacement_fluct(size,F) - outer = np.dot(np.average(F,axis=(0,1,2)),size) + 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] + 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) + tree = _spatial.cKDTree(c.reshape(-1,3),boxsize=outer) return tree.query(cell_coord0(new_grid,outer))[1].flatten() diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 674ff9c5a..e19f140fb 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,4 +1,4 @@ -import numpy as np +import numpy as _np def Cauchy(P,F): """ @@ -14,10 +14,10 @@ def Cauchy(P,F): First Piola-Kirchhoff stress. """ - if np.shape(F) == np.shape(P) == (3,3): - sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) + if _np.shape(F) == _np.shape(P) == (3,3): + sigma = 1.0/_np.linalg.det(F) * _np.dot(P,F.T) else: - sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F) + sigma = _np.einsum('i,ijk,ilk->ijl',1.0/_np.linalg.det(F),P,F) return symmetric(sigma) @@ -31,8 +31,8 @@ def deviatoric_part(T): Tensor of which the deviatoric part is computed. """ - return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \ - T - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[T.shape[0],3,3]),spherical_part(T)) + return T - _np.eye(3)*spherical_part(T) if _np.shape(T) == (3,3) else \ + T - _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),[T.shape[0],3,3]),spherical_part(T)) def eigenvalues(T_sym): @@ -48,7 +48,7 @@ def eigenvalues(T_sym): Symmetric tensor of which the eigenvalues are computed. """ - return np.linalg.eigvalsh(symmetric(T_sym)) + return _np.linalg.eigvalsh(symmetric(T_sym)) def eigenvectors(T_sym,RHS=False): @@ -65,13 +65,13 @@ def eigenvectors(T_sym,RHS=False): Enforce right-handed coordinate system. Default is False. """ - (u,v) = np.linalg.eigh(symmetric(T_sym)) + (u,v) = _np.linalg.eigh(symmetric(T_sym)) if RHS: - if np.shape(T_sym) == (3,3): - if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 + if _np.shape(T_sym) == (3,3): + if _np.linalg.det(v) < 0.0: v[:,2] *= -1.0 else: - v[np.linalg.det(v) < 0.0,:,2] *= -1.0 + v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 return v @@ -99,7 +99,7 @@ def maximum_shear(T_sym): """ w = eigenvalues(T_sym) - return (w[0] - w[2])*0.5 if np.shape(T_sym) == (3,3) else \ + return (w[0] - w[2])*0.5 if _np.shape(T_sym) == (3,3) else \ (w[:,0] - w[:,2])*0.5 @@ -141,10 +141,10 @@ def PK2(P,F): Deformation gradient. """ - if np.shape(F) == np.shape(P) == (3,3): - S = np.dot(np.linalg.inv(F),P) + if _np.shape(F) == _np.shape(P) == (3,3): + S = _np.dot(_np.linalg.inv(F),P) else: - S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P) + S = _np.einsum('ijk,ikl->ijl',_np.linalg.inv(F),P) return symmetric(S) @@ -187,14 +187,14 @@ def spherical_part(T,tensor=False): """ if T.shape == (3,3): - sph = np.trace(T)/3.0 - return sph if not tensor else np.eye(3)*sph + sph = _np.trace(T)/3.0 + return sph if not tensor else _np.eye(3)*sph else: - sph = np.trace(T,axis1=1,axis2=2)/3.0 + sph = _np.trace(T,axis1=1,axis2=2)/3.0 if not tensor: return sph else: - return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(T.shape[0],3,3)),sph) + return _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),(T.shape[0],3,3)),sph) def strain_tensor(F,t,m): @@ -216,22 +216,22 @@ def strain_tensor(F,t,m): """ F_ = F.reshape(1,3,3) if F.shape == (3,3) else F if t == 'V': - B = np.matmul(F_,transpose(F_)) - w,n = np.linalg.eigh(B) + B = _np.matmul(F_,transpose(F_)) + w,n = _np.linalg.eigh(B) elif t == 'U': - C = np.matmul(transpose(F_),F_) - w,n = np.linalg.eigh(C) + C = _np.matmul(transpose(F_),F_) + w,n = _np.linalg.eigh(C) if m > 0.0: - eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - - np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) + - _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3])) elif m < 0.0: - eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - + np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) + + _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3])) else: - eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) + eps = _np.matmul(n,_np.einsum('ij,ikj->ijk',0.5*_np.log(w),n)) - return eps.reshape(3,3) if np.shape(F) == (3,3) else \ + return eps.reshape(3,3) if _np.shape(F) == (3,3) else \ eps @@ -258,8 +258,8 @@ def transpose(T): Tensor of which the transpose is computed. """ - return T.T if np.shape(T) == (3,3) else \ - np.transpose(T,(0,2,1)) + return T.T if _np.shape(T) == (3,3) else \ + _np.transpose(T,(0,2,1)) def _polar_decomposition(T,requested): @@ -275,17 +275,17 @@ def _polar_decomposition(T,requested): ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. """ - u, s, vh = np.linalg.svd(T) - R = np.dot(u,vh) if np.shape(T) == (3,3) else \ - np.einsum('ijk,ikl->ijl',u,vh) + u, s, vh = _np.linalg.svd(T) + R = _np.dot(u,vh) if _np.shape(T) == (3,3) else \ + _np.einsum('ijk,ikl->ijl',u,vh) output = [] if 'R' in requested: output.append(R) if 'V' in requested: - output.append(np.dot(T,R.T) if np.shape(T) == (3,3) else np.einsum('ijk,ilk->ijl',T,R)) + output.append(_np.dot(T,R.T) if _np.shape(T) == (3,3) else _np.einsum('ijk,ilk->ijl',T,R)) if 'U' in requested: - output.append(np.dot(R.T,T) if np.shape(T) == (3,3) else np.einsum('ikj,ikl->ijl',R,T)) + output.append(_np.dot(R.T,T) if _np.shape(T) == (3,3) else _np.einsum('ikj,ikl->ijl',R,T)) return tuple(output) @@ -303,5 +303,5 @@ def _Mises(T_sym,s): """ d = deviatoric_part(T_sym) - return np.sqrt(s*(np.sum(d**2.0))) if np.shape(T_sym) == (3,3) else \ - np.sqrt(s*np.einsum('ijk->i',d**2.0)) + return _np.sqrt(s*(_np.sum(d**2.0))) if _np.shape(T_sym) == (3,3) else \ + _np.sqrt(s*_np.einsum('ijk->i',d**2.0)) diff --git a/python/damask/util.py b/python/damask/util.py index 273da8c1e..d45ea366e 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -9,37 +9,22 @@ from optparse import Option import numpy as np -class bcolors: - """ - ASCII Colors. - - https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py - https://stackoverflow.com/questions/287871 - """ - - HEADER = '\033[95m' - OKBLUE = '\033[94m' - OKGREEN = '\033[92m' - WARNING = '\033[93m' - FAIL = '\033[91m' - ENDC = '\033[0m' - BOLD = '\033[1m' - DIM = '\033[2m' - UNDERLINE = '\033[4m' - CROSSOUT = '\033[9m' - - def disable(self): - self.HEADER = '' - self.OKBLUE = '' - self.OKGREEN = '' - self.WARNING = '' - self.FAIL = '' - self.ENDC = '' - self.BOLD = '' - self.UNDERLINE = '' - self.CROSSOUT = '' - +# limit visibility +__all__=[ + 'srepr', + 'croak', + 'report', + 'emph','deemph','delete','strikeout', + 'execute', + 'show_progress', + 'scale_to_coprime', + 'return_message', + 'extendableOption', + ] +#################################################################################################### +# Functions +#################################################################################################### def srepr(arg,glue = '\n'): r""" Join arguments as individual lines. @@ -144,6 +129,52 @@ def execute(cmd, return out,error +def show_progress(iterable,N_iter=None,prefix='',bar_length=50): + """ + Decorate a loop with a status bar. + + Use similar like enumerate. + + Parameters + ---------- + iterable : iterable/function with yield statement + Iterable (or function with yield statement) to be decorated. + N_iter : int + Total # of iterations. Needed if number of iterations can not be obtained as len(iterable). + prefix : str, optional. + Prefix string. + bar_length : int, optional + Character length of bar. Defaults to 50. + + """ + status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length) + + for i,item in enumerate(iterable): + yield item + status.update(i) + + +def scale_to_coprime(v): + """Scale vector to co-prime (relatively prime) integers.""" + MAX_DENOMINATOR = 1000 + + def get_square_denominator(x): + """Denominator of the square of a number.""" + return fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator + + def lcm(a, b): + """Least common multiple.""" + return a * b // np.gcd(a, b) + + denominators = [int(get_square_denominator(i)) for i in v] + s = reduce(lcm, denominators) ** 0.5 + m = (np.array(v)*s).astype(np.int) + return m//reduce(np.gcd,m) + + +#################################################################################################### +# Classes +#################################################################################################### class extendableOption(Option): """ Used for definition of new option parser action 'extend', which enables to take multiple option arguments. @@ -215,47 +246,36 @@ class _ProgressBar: sys.stderr.write('\n') sys.stderr.flush() -def show_progress(iterable,N_iter=None,prefix='',bar_length=50): + +class bcolors: """ - Decorate a loop with a status bar. - - Use similar like enumerate. - - Parameters - ---------- - iterable : iterable/function with yield statement - Iterable (or function with yield statement) to be decorated. - N_iter : int - Total # of iterations. Needed if number of iterations can not be obtained as len(iterable). - prefix : str, optional. - Prefix string. - bar_length : int, optional - Character length of bar. Defaults to 50. + ASCII Colors. + https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py + https://stackoverflow.com/questions/287871 """ - status = _ProgressBar(N_iter if N_iter else len(iterable),prefix,bar_length) - for i,item in enumerate(iterable): - yield item - status.update(i) + HEADER = '\033[95m' + OKBLUE = '\033[94m' + OKGREEN = '\033[92m' + WARNING = '\033[93m' + FAIL = '\033[91m' + ENDC = '\033[0m' + BOLD = '\033[1m' + DIM = '\033[2m' + UNDERLINE = '\033[4m' + CROSSOUT = '\033[9m' - -def scale_to_coprime(v): - """Scale vector to co-prime (relatively prime) integers.""" - MAX_DENOMINATOR = 1000 - - def get_square_denominator(x): - """Denominator of the square of a number.""" - return fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator - - def lcm(a, b): - """Least common multiple.""" - return a * b // np.gcd(a, b) - - denominators = [int(get_square_denominator(i)) for i in v] - s = reduce(lcm, denominators) ** 0.5 - m = (np.array(v)*s).astype(np.int) - return m//reduce(np.gcd,m) + def disable(self): + self.HEADER = '' + self.OKBLUE = '' + self.OKGREEN = '' + self.WARNING = '' + self.FAIL = '' + self.ENDC = '' + self.BOLD = '' + self.UNDERLINE = '' + self.CROSSOUT = '' class return_message: @@ -276,4 +296,3 @@ class return_message: def __repr__(self): """Return message suitable for interactive shells.""" return srepr(self.message) - diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index cdddca89e..acbdbf688 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -24,7 +24,7 @@ class TestGridFilters: n = grid_filters.node_coord0(grid,size) + size/grid*.5 assert np.allclose(c,n) - @pytest.mark.parametrize('mode',[('cell'),('node')]) + @pytest.mark.parametrize('mode',['cell','node']) def test_grid_DNA(self,mode): """Ensure that xx_coord0_gridSizeOrigin is the inverse of xx_coord0.""" grid = np.random.randint(8,32,(3)) @@ -49,7 +49,7 @@ class TestGridFilters: assert np.allclose(grid_filters.node_coord(size,F) [1:-1,1:-1,1:-1],grid_filters.cell_2_node( grid_filters.cell_coord(size,F))[1:-1,1:-1,1:-1]) - @pytest.mark.parametrize('mode',[('cell'),('node')]) + @pytest.mark.parametrize('mode',['cell','node']) def test_coord0_origin(self,mode): origin= np.random.random(3) size = np.random.random(3) # noqa @@ -61,22 +61,24 @@ class TestGridFilters: elif mode == 'node': assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,))) - @pytest.mark.parametrize('mode',[('cell'),('node')]) - def test_displacement_avg_vanishes(self,mode): + @pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg, + grid_filters.node_displacement_avg]) + def test_displacement_avg_vanishes(self,function): """Ensure that random fluctuations in F do not result in average displacement.""" - size = np.random.random(3) # noqa + 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) + assert np.allclose(function(size,F),0.0) - @pytest.mark.parametrize('mode',[('cell'),('node')]) - def test_displacement_fluct_vanishes(self,mode): + @pytest.mark.parametrize('function',[grid_filters.cell_displacement_fluct, + grid_filters.node_displacement_fluct]) + def test_displacement_fluct_vanishes(self,function): """Ensure that constant F does not result in fluctuating displacement.""" - size = np.random.random(3) # noqa + 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)) # noqa - assert np.allclose(eval('grid_filters.{}_displacement_fluct(size,F)'.format(mode)),0.0) + F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) + assert np.allclose(function(size,F),0.0) def test_regrid(self): size = np.random.random(3) diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90 new file mode 100644 index 000000000..7d3043ed0 --- /dev/null +++ b/src/LAPACK_interface.f90 @@ -0,0 +1,59 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Fortran interfaces for LAPACK routines +!> @details https://www.netlib.org/lapack/ +!-------------------------------------------------------------------------------------------------- +module LAPACK_interface + interface + + subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info) + use prec + character, intent(in) :: jobvl,jobvr + integer, intent(in) :: n,lda,ldvl,ldvr,lwork + real(pReal), intent(inout), dimension(lda,n) :: a + real(pReal), intent(out), dimension(n) :: wr,wi + real(pReal), intent(out), dimension(ldvl,n) :: vl + real(pReal), intent(out), dimension(ldvr,n) :: vr + real(pReal), intent(out), dimension(max(1,lwork)) :: work + integer, intent(out) :: info + end subroutine dgeev + + subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info) + use prec + integer, intent(in) :: n,nrhs,lda,ldb + real(pReal), intent(inout), dimension(lda,n) :: a + integer, intent(out), dimension(n) :: ipiv + real(pReal), intent(out), dimension(ldb,nrhs) :: b + integer, intent(out) :: info + end subroutine dgesv + + subroutine dgetrf(m,n,a,lda,ipiv,info) + use prec + integer, intent(in) :: m,n,lda + real(pReal), intent(inout), dimension(lda,n) :: a + integer, intent(out), dimension(min(m,n)) :: ipiv + integer, intent(out) :: info + end subroutine dgetrf + + subroutine dgetri(n,a,lda,ipiv,work,lwork,info) + use prec + integer, intent(in) :: n,lda,lwork + real(pReal), intent(inout), dimension(lda,n) :: a + integer, intent(out), dimension(n) :: ipiv + real(pReal), intent(out), dimension(max(1,lwork)) :: work + integer, intent(out) :: info + end subroutine dgetri + + subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info) + use prec + character, intent(in) :: jobz,uplo + integer, intent(in) :: n,lda,lwork + real(pReal), intent(inout), dimension(lda,n) :: a + real(pReal), intent(out), dimension(n) :: w + real(pReal), intent(out), dimension(max(1,lwork)) :: work + integer, intent(out) :: info + end subroutine dsyev + + end interface + +end module LAPACK_interface diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index ab26ee9d5..64ad3e1d7 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -9,6 +9,7 @@ #include "list.f90" #include "future.f90" #include "config.f90" +#include "LAPACK_interface.f90" #include "math.f90" #include "quaternions.f90" #include "Lambert.f90" diff --git a/src/crystallite.f90 b/src/crystallite.f90 index efa066696..9bc254e0c 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -835,8 +835,6 @@ logical function integrateStress(ipc,ip,el,timeFraction) jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update logical :: error - external :: & - dgesv integrateStress = .false. diff --git a/src/math.f90 b/src/math.f90 index b0852e8d4..070751de8 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -9,6 +9,7 @@ module math use prec use IO use numerics + use LAPACK_interface implicit none public @@ -485,18 +486,14 @@ function math_invSym3333(A) real(pReal),dimension(3,3,3,3),intent(in) :: A - integer, dimension(6) :: ipiv6 - real(pReal), dimension(6,6) :: temp66 - real(pReal), dimension(6*(64+2)) :: work - integer :: ierr_i, ierr_f - external :: & - dgetrf, & - dgetri + integer, dimension(6) :: ipiv6 + real(pReal), dimension(6,6) :: temp66 + real(pReal), dimension(6*6) :: work + integer :: ierr_i, ierr_f temp66 = math_sym3333to66(A) call dgetrf(6,6,temp66,6,ipiv6,ierr_i) call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f) - if (ierr_i /= 0 .or. ierr_f /= 0) then call IO_error(400, ext_msg = 'math_invSym3333') else @@ -515,12 +512,9 @@ subroutine math_invert(InvA, error, A) real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA logical, intent(out) :: error - integer, dimension(size(A,1)) :: ipiv - real(pReal), dimension(size(A,1)*(64+2)) :: work - integer :: ierr - external :: & - dgetrf, & - dgetri + integer, dimension(size(A,1)) :: ipiv + real(pReal), dimension(size(A,1)**2) :: work + integer :: ierr invA = A call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr) @@ -884,9 +878,7 @@ subroutine math_eigh(m,w,v,error) logical, intent(out) :: error integer :: ierr - real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f - external :: & - dsyev + real(pReal), dimension(size(m,1)**2) :: work v = m ! copy matrix to input (doubles as output) array call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr) @@ -1041,9 +1033,7 @@ function math_eigvalsh(m) real(pReal), dimension(size(m,1),size(m,1)) :: m_ integer :: ierr - real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f - external :: & - dsyev + real(pReal), dimension(size(m,1)**2) :: work m_= m ! copy matrix to input (will be destroyed) call dsyev('N','U',size(m,1),m_,size(m,1),math_eigvalsh,work,size(work,1),ierr) diff --git a/src/rotations.f90 b/src/rotations.f90 index 7ce366f74..b4e143f28 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -596,8 +596,6 @@ function om2ax(om) result(ax) real(pReal), dimension(3,3) :: VR, devNull, om_ integer :: ierr, i - external :: dgeev - om_ = om ! first get the rotation angle