223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'vector-mechanics' into 'development' 

Vector mechanics

See merge request damask/DAMASK!281
This commit is contained in:
Sharan Roongta 2020-11-24 11:41:45 +01:00
parent 9ef19f5b1a 5ab2847b36
commit 6e1906eae0
19 changed files with 748 additions and 465 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 25b7b78951ab9fb682aaf048d3e0f0d09a3be695
Subproject commit 8bd09b5511d1e0e0ea288b47d16ce4924d75adcd

14
processing/post/addOrientations.py
View File

@ -115,10 +115,10 @@ for name in filenames:
if options.eulers is not None:
label = options.eulers
print(np.max(table.get(options.eulers),axis=0))
o = damask.Rotation.from_Eulers(table.get(options.eulers), options.degrees)
o = damask.Rotation.from_Euler_angles(table.get(options.eulers), options.degrees)
elif options.rodrigues is not None:
label = options.rodrigues
o = damask.Rotation.from_Rodrigues(table.get(options.rodrigues))
o = damask.Rotation.from_Rodrigues_vector(table.get(options.rodrigues))
elif options.matrix is not None:
label = options.matrix
o = damask.Rotation.from_matrix(table.get(options.matrix).reshape(-1,3,3))
@ -137,14 +137,14 @@ for name in filenames:
if 'rodrigues' in options.output:
table = table.add('ro({})'.format(label),o.as_Rodrigues(), scriptID+' '+' '.join(sys.argv[1:]))
table = table.add('ro({})'.format(label),o.as_Rodrigues_vector(), scriptID+' '+' '.join(sys.argv[1:]))
if 'eulers' in options.output:
table = table.add('eu({})'.format(label),o.as_Eulers(options.degrees), scriptID+' '+' '.join(sys.argv[1:]))
table = table.add('eu({})'.format(label),o.as_Euler_angles(options.degrees),scriptID+' '+' '.join(sys.argv[1:]))
if 'quaternion' in options.output:
table = table.add('qu({})'.format(label),o.as_quaternion(), scriptID+' '+' '.join(sys.argv[1:]))
table = table.add('qu({})'.format(label),o.as_quaternion(), scriptID+' '+' '.join(sys.argv[1:]))
if 'matrix' in options.output:
table = table.add('om({})'.format(label),o.as_matrix(), scriptID+' '+' '.join(sys.argv[1:]))
table = table.add('om({})'.format(label),o.as_matrix(), scriptID+' '+' '.join(sys.argv[1:]))
if 'axisangle' in options.output:
table = table.add('om({})'.format(label),o.as_axisangle(options.degrees), scriptID+' '+' '.join(sys.argv[1:]))
table = table.add('om({})'.format(label),o.as_axis_angle(options.degrees), scriptID+' '+' '.join(sys.argv[1:]))
table.save((sys.stdout if name is None else name), legacy=True)

11
python/damask/__init__.py
View File

@ -1,4 +1,12 @@
"""Tools for pre and post processing of DAMASK simulations."""
"""
Tools for pre and post processing of DAMASK simulations.
Modules that contain only one class (of the same name),
are prefixed by a '_'. For example, '_colormap' contains
a class called 'Colormap' which is imported as 'damask.Colormap'.
"""
from pathlib import Path as _Path
import re as _re
@ -12,6 +20,7 @@ from ._environment import Environment as _ # noqa
environment = _()
from . import util # noqa
from . import seeds # noqa
from . import tensor # noqa
from . import mechanics # noqa
from . import solver # noqa
from . import grid_filters # noqa

2
python/damask/_geom.py
View File

@ -759,7 +759,7 @@ class Geom:
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
Eulers = R.as_Eulers(degrees=True)
Eulers = R.as_Euler_angles(degrees=True)
material_in = self.material.copy()
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')

40
python/damask/_orientation.py
View File

@ -2,7 +2,7 @@ import numpy as np
from . import Rotation
from . import util
from . import mechanics
from . import tensor
_parameter_doc = \
"""lattice : str
@ -67,6 +67,7 @@ class Orientation(Rotation):
Examples
--------
An array of 3 x 5 random orientations reduced to the fundamental zone of tetragonal symmetry:
>>> damask.Orientation.from_random(shape=(3,5),lattice='tetragonal').reduced
"""
@ -263,9 +264,9 @@ class Orientation(Rotation):
@classmethod
@util.extended_docstring(Rotation.from_Eulers,_parameter_doc)
def from_Eulers(cls,**kwargs):
return cls(rotation=Rotation.from_Eulers(**kwargs),**kwargs)
@util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
def from_Euler_angles(cls,**kwargs):
return cls(rotation=Rotation.from_Euler_angles(**kwargs),**kwargs)
@classmethod
@ -287,9 +288,9 @@ class Orientation(Rotation):
@classmethod
@util.extended_docstring(Rotation.from_Rodrigues,_parameter_doc)
def from_Rodrigues(cls,**kwargs):
return cls(rotation=Rotation.from_Rodrigues(**kwargs),**kwargs)
@util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
def from_Rodrigues_vector(cls,**kwargs):
return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs),**kwargs)
@classmethod
@ -334,7 +335,7 @@ class Orientation(Rotation):
x = o.to_frame(uvw=uvw)
z = o.to_frame(hkl=hkl)
om = np.stack([x,np.cross(z,x),z],axis=-2)
return o.copy(rotation=Rotation.from_matrix(mechanics.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
@property
@ -468,7 +469,7 @@ class Orientation(Rotation):
if self.family is None:
raise ValueError('Missing crystal symmetry')
rho_abs = np.abs(self.as_Rodrigues(vector=True))
rho_abs = np.abs(self.as_Rodrigues_vector(compact=True))
with np.errstate(invalid='ignore'):
# using '*'/prod for 'and'
@ -511,7 +512,7 @@ class Orientation(Rotation):
if self.family is None:
raise ValueError('Missing crystal symmetry')
rho = self.as_Rodrigues(vector=True)
rho = self.as_Rodrigues_vector(compact=True)
with np.errstate(invalid='ignore'):
if self.family == 'cubic':
@ -856,7 +857,6 @@ class Orientation(Rotation):
vector_ = self.to_SST(vector,proper) if in_SST else \
self @ np.broadcast_to(vector,self.shape+(3,))
if self.family == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
@ -962,9 +962,9 @@ class Orientation(Rotation):
"""
if self.family is None or other.family is None:
raise ValueError('Missing crystal symmetry')
raise ValueError('missing crystal symmetry')
if self.family != other.family:
raise NotImplementedError('Disorientation between different crystal families not supported yet.')
raise NotImplementedError('disorientation between different crystal families')
blend = util.shapeblender(self.shape,other.shape)
s = self.equivalent
@ -1106,7 +1106,7 @@ class Orientation(Rotation):
(np.array(hkil),np.array([[1,0,0,0],
[0,1,0,0],
[0,0,0,1]]))
return np.einsum('il,...l->...i',basis,axis)
return np.einsum('il,...l',basis,axis)
@classmethod
@ -1136,7 +1136,7 @@ class Orientation(Rotation):
[ 0, 1, 0],
[-1,-1, 0],
[ 0, 0, 1]]))
return np.einsum('il,...l->...i',basis,axis)
return np.einsum('il,...l',basis,axis)
def to_lattice(self,*,direction=None,plane=None):
@ -1160,7 +1160,7 @@ class Orientation(Rotation):
axis,basis = (np.array(direction),self.basis_reciprocal.T) \
if plane is None else \
(np.array(plane),self.basis_real.T)
return np.einsum('il,...l->...i',basis,axis)
return np.einsum('il,...l',basis,axis)
def to_frame(self,*,uvw=None,hkl=None,with_symmetry=False):
@ -1186,9 +1186,9 @@ class Orientation(Rotation):
if hkl is None else \
(np.array(hkl),self.basis_reciprocal)
return (self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+axis.shape[:-1],mode='right')
@ np.broadcast_to(np.einsum('il,...l->...i',basis,axis),self.symmetry_operations.shape+axis.shape)
@ np.broadcast_to(np.einsum('il,...l',basis,axis),self.symmetry_operations.shape+axis.shape)
if with_symmetry else
np.einsum('il,...l->...i',basis,axis))
np.einsum('il,...l',basis,axis))
def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False):
@ -1243,8 +1243,8 @@ class Orientation(Rotation):
"""
d = self.to_frame(uvw=self.kinematics[mode]['direction'],with_symmetry=False)
p = self.to_frame(hkl=self.kinematics[mode]['plane'] ,with_symmetry=False)
P = np.einsum('...i,...j->...ij',d/np.linalg.norm(d,axis=-1,keepdims=True),
p/np.linalg.norm(p,axis=-1,keepdims=True))
P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=-1,keepdims=True),
p/np.linalg.norm(p,axis=-1,keepdims=True))
return ~self.broadcast_to( self.shape+P.shape[:-2],mode='right') \
@ np.broadcast_to(P,self.shape+P.shape)

62
python/damask/_result.py
View File

@ -18,6 +18,7 @@ from . import Table
from . import Orientation
from . import grid_filters
from . import mechanics
from . import tensor
from . import util
h5py3 = h5py.__version__[0] == '3'
@ -588,19 +589,19 @@ class Result:
@staticmethod
def _add_Cauchy(P,F):
def _add_stress_Cauchy(P,F):
return {
'data': mechanics.Cauchy(P['data'],F['data']),
'data': mechanics.stress_Cauchy(P['data'],F['data']),
'label': 'sigma',
'meta': {
'Unit': P['meta']['Unit'],
'Description': "Cauchy stress calculated "
f"from {P['label']} ({P['meta']['Description']})"
f" and {F['label']} ({F['meta']['Description']})",
'Creator': 'add_Cauchy'
'Creator': 'add_stress_Cauchy'
}
}
def add_Cauchy(self,P='P',F='F'):
def add_stress_Cauchy(self,P='P',F='F'):
"""
Add Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
@ -612,7 +613,7 @@ class Result:
Label of the dataset containing the deformation gradient. Defaults to ‘F’.
"""
self._add_generic_pointwise(self._add_Cauchy,{'P':P,'F':F})
self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F})
@staticmethod
@ -642,7 +643,7 @@ class Result:
@staticmethod
def _add_deviator(T):
return {
'data': mechanics.deviatoric_part(T['data']),
'data': tensor.deviatoric(T['data']),
'label': f"s_{T['label']}",
'meta': {
'Unit': T['meta']['Unit'],
@ -673,7 +674,7 @@ class Result:
label,p = 'Minimum',0
return {
'data': mechanics.eigenvalues(T_sym['data'])[:,p],
'data': tensor.eigenvalues(T_sym['data'])[:,p],
'label': f"lambda_{eigenvalue}({T_sym['label']})",
'meta' : {
'Unit': T_sym['meta']['Unit'],
@ -705,7 +706,7 @@ class Result:
elif eigenvalue == 'min':
label,p = 'minimum',0
return {
'data': mechanics.eigenvectors(T_sym['data'])[:,p],
'data': tensor.eigenvectors(T_sym['data'])[:,p],
'label': f"v_{eigenvalue}({T_sym['label']})",
'meta' : {
'Unit': '1',
@ -789,7 +790,7 @@ class Result:
@staticmethod
def _add_Mises(T_sym,kind):
def _add_equivalent_Mises(T_sym,kind):
k = kind
if k is None:
if T_sym['meta']['Unit'] == '1':
@ -800,7 +801,8 @@ class Result:
raise ValueError('invalid von Mises kind {kind}')
return {
'data': (mechanics.Mises_strain if k=='strain' else mechanics.Mises_stress)(T_sym['data']),
'data': (mechanics.equivalent_strain_Mises if k=='strain' else \
mechanics.equivalent_stress_Mises)(T_sym['data']),
'label': f"{T_sym['label']}_vM",
'meta': {
'Unit': T_sym['meta']['Unit'],
@ -808,7 +810,7 @@ class Result:
'Creator': 'add_Mises'
}
}
def add_Mises(self,T_sym,kind=None):
def add_equivalent_Mises(self,T_sym,kind=None):
"""
Add the equivalent Mises stress or strain of a symmetric tensor.
@ -821,7 +823,7 @@ class Result:
it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress').
"""
self._add_generic_pointwise(self._add_Mises,{'T_sym':T_sym},{'kind':kind})
self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind})
@staticmethod
@ -863,19 +865,19 @@ class Result:
@staticmethod
def _add_PK2(P,F):
def _add_stress_second_Piola_Kirchhoff(P,F):
return {
'data': mechanics.PK2(P['data'],F['data']),
'data': mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']),
'label': 'S',
'meta': {
'Unit': P['meta']['Unit'],
'Description': "2. Piola-Kirchhoff stress calculated "
f"from {P['label']} ({P['meta']['Description']})"
f" and {F['label']} ({F['meta']['Description']})",
'Creator': 'add_PK2'
'Creator': 'add_stress_second_Piola_Kirchhoff'
}
}
def add_PK2(self,P='P',F='F'):
def add_stress_second_Piola_Kirchhoff(self,P='P',F='F'):
"""
Add second Piola-Kirchhoff stress calculated from first Piola-Kirchhoff stress and deformation gradient.
@ -887,7 +889,7 @@ class Result:
Label of deformation gradient dataset. Defaults to ‘F’.
"""
self._add_generic_pointwise(self._add_PK2,{'P':P,'F':F})
self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F})
# The add_pole functionality needs discussion.
@ -934,17 +936,17 @@ class Result:
@staticmethod
def _add_rotational_part(F):
def _add_rotation(F):
return {
'data': mechanics.rotational_part(F['data']),
'data': mechanics.rotation(F['data']).as_matrix(),
'label': f"R({F['label']})",
'meta': {
'Unit': F['meta']['Unit'],
'Description': f"Rotational part of {F['label']} ({F['meta']['Description']})",
'Creator': 'add_rotational_part'
'Creator': 'add_rotation'
}
}
def add_rotational_part(self,F):
def add_rotation(self,F):
"""
Add rotational part of a deformation gradient.
@ -954,13 +956,13 @@ class Result:
Label of deformation gradient dataset.
"""
self._add_generic_pointwise(self._add_rotational_part,{'F':F})
self._add_generic_pointwise(self._add_rotation,{'F':F})
@staticmethod
def _add_spherical(T):
return {
'data': mechanics.spherical_part(T['data']),
'data': tensor.spherical(T['data'],False),
'label': f"p_{T['label']}",
'meta': {
'Unit': T['meta']['Unit'],
@ -982,21 +984,21 @@ class Result:
@staticmethod
def _add_strain_tensor(F,t,m):
def _add_strain(F,t,m):
return {
'data': mechanics.strain_tensor(F['data'],t,m),
'data': mechanics.strain(F['data'],t,m),
'label': f"epsilon_{t}^{m}({F['label']})",
'meta': {
'Unit': F['meta']['Unit'],
'Description': f"Strain tensor of {F['label']} ({F['meta']['Description']})",
'Creator': 'add_strain_tensor'
'Creator': 'add_strain'
}
}
def add_strain_tensor(self,F='F',t='V',m=0.0):
def add_strain(self,F='F',t='V',m=0.0):
"""
Add strain tensor of a deformation gradient.
For details refer to damask.mechanics.strain_tensor
For details refer to damask.mechanics.strain
Parameters
----------
@ -1009,13 +1011,13 @@ class Result:
Order of the strain calculation. Defaults to ‘0.0’.
"""
self._add_generic_pointwise(self._add_strain_tensor,{'F':F},{'t':t,'m':m})
self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m})
@staticmethod
def _add_stretch_tensor(F,t):
return {
'data': (mechanics.left_stretch if t.upper() == 'V' else mechanics.right_stretch)(F['data']),
'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']),
'label': f"{t}({F['label']})",
'meta': {
'Unit': F['meta']['Unit'],

66
python/damask/_rotation.py
View File

@ -1,6 +1,6 @@
import numpy as np
from . import mechanics
from . import tensor
from . import util
from . import grid_filters
@ -61,18 +61,21 @@ class Rotation:
elif np.array(rotation).shape[-1] == 4:
self.quaternion = np.array(rotation)
else:
raise ValueError('"rotation" is neither a Rotation nor a quaternion')
raise TypeError('"rotation" is neither a Rotation nor a quaternion')
def __repr__(self):
"""Represent rotation as unit quaternion, rotation matrix, and Bunge-Euler angles."""
return 'Quaternions:\n'+str(self.quaternion) \
if self.quaternion.shape != (4,) else \
'\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)),
])
if self == Rotation():
return 'Rotation()'
else:
return f'Quaternions {self.shape}:\n'+str(self.quaternion) \
if self.quaternion.shape != (4,) else \
'\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Euler_angles(degrees=True)),
])
# ToDo: Check difference __copy__ vs __deepcopy__
@ -298,8 +301,8 @@ class Rotation:
"""
return self.quaternion.copy()
def as_Eulers(self,
degrees = False):
def as_Euler_angles(self,
degrees = False):
"""
Represent as Bunge-Euler angles.
@ -356,8 +359,8 @@ class Rotation:
"""
return Rotation._qu2om(self.quaternion)
def as_Rodrigues(self,
vector = False):
def as_Rodrigues_vector(self,
compact = False):
"""
Represent as Rodrigues-Frank vector with separated axis and angle argument.
@ -375,7 +378,7 @@ class Rotation:
"""
ro = Rotation._qu2ro(self.quaternion)
if vector:
if compact:
with np.errstate(invalid='ignore'):
return ro[...,:3]*ro[...,3:4]
else:
@ -446,9 +449,9 @@ class Rotation:
return Rotation(qu)
@staticmethod
def from_Eulers(phi,
degrees = False,
**kwargs):
def from_Euler_angles(phi,
degrees = False,
**kwargs):
"""
Initialize from Bunge-Euler angles.
@ -533,11 +536,11 @@ class Rotation:
raise ValueError('Invalid shape.')
if reciprocal:
om = np.linalg.inv(mechanics.transpose(om)/np.pi) # transform reciprocal basis set
om = np.linalg.inv(tensor.transpose(om)/np.pi) # transform reciprocal basis set
orthonormal = False # contains stretch
if not orthonormal:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.einsum('...ij,...jl->...il',U,Vh)
om = np.einsum('...ij,...jl',U,Vh)
if not np.all(np.isclose(np.linalg.det(om),1.0)):
raise ValueError('Orientation matrix has determinant ≠ 1.')
if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \
@ -590,18 +593,17 @@ class Rotation:
@staticmethod
def from_Rodrigues(rho,
normalize = False,
P = -1,
**kwargs):
def from_Rodrigues_vector(rho,
normalize = False,
P = -1,
**kwargs):
"""
Initialize from Rodrigues-Frank vector.
Initialize from Rodrigues-Frank vector (angle separated from axis).
Parameters
----------
rho : numpy.ndarray of shape (...,4)
Rodrigues-Frank vector (angle separated from axis).
(n_1, n_2, n_3, tan(ω/2)), ǀnǀ = 1 and ω ∈ [0,π].
Rodrigues-Frank vector. (n_1, n_2, n_3, tan(ω/2)), ǀnǀ = 1 and ω ∈ [0,π].
normalize : boolean, optional
Allow ǀnǀ ≠ 1. Defaults to False.
P : int ∈ {-1,1}, optional
@ -714,7 +716,7 @@ class Rotation:
@staticmethod
def from_ODF(weights,
Eulers,
phi,
N = 500,
degrees = True,
fractions = True,
@ -727,7 +729,7 @@ class Rotation:
----------
weights : numpy.ndarray of shape (n)
Texture intensity values (probability density or volume fraction) at Euler grid points.
Eulers : numpy.ndarray of shape (n,3)
phi : numpy.ndarray of shape (n,3)
Grid coordinates in Euler space at which weights are defined.
N : integer, optional
Number of discrete orientations to be sampled from the given ODF.
@ -760,15 +762,15 @@ class Rotation:
"""
def _dg(eu,deg):
"""Return infinitesimal Euler space volume of bin(s)."""
Eulers_sorted = eu[np.lexsort((eu[:,0],eu[:,1],eu[:,2]))]
steps,size,_ = grid_filters.cell_coord0_gridSizeOrigin(Eulers_sorted)
phi_sorted = eu[np.lexsort((eu[:,0],eu[:,1],eu[:,2]))]
steps,size,_ = grid_filters.cell_coord0_gridSizeOrigin(phi_sorted)
delta = np.radians(size/steps) if deg else size/steps
return delta[0]*2.0*np.sin(delta[1]/2.0)*delta[2] / 8.0 / np.pi**2 * np.sin(np.radians(eu[:,1]) if deg else eu[:,1])
dg = 1.0 if fractions else _dg(Eulers,degrees)
dg = 1.0 if fractions else _dg(phi,degrees)
dV_V = dg * np.maximum(0.0,weights.squeeze())
return Rotation.from_Eulers(Eulers[util.hybrid_IA(dV_V,N,rng_seed)],degrees)
return Rotation.from_Euler_angles(phi[util.hybrid_IA(dV_V,N,rng_seed)],degrees)
@staticmethod

3
python/damask/_vtk.py
View File

@ -1,3 +1,4 @@
import os
import multiprocessing as mp
from pathlib import Path
@ -138,6 +139,8 @@ class VTK:
vtkUnstructuredGrid, and vtkPolyData.
"""
if not os.path.isfile(fname): # vtk has a strange error handling
raise FileNotFoundError(f'no such file: {fname}')
ext = Path(fname).suffix
if ext == '.vtk' or dataset_type is not None:
reader = vtk.vtkGenericDataObjectReader()

343
python/damask/mechanics.py
View File

@ -1,307 +1,282 @@
"""Finite-strain continuum mechanics."""
from . import tensor as _tensor
from . import _rotation
import numpy as _np
def Cauchy(P,F):
"""
Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
def deformation_Cauchy_Green_left(F):
"""
Calculate left Cauchy-Green deformation tensor (Finger deformation tensor).
Parameters
----------
F : numpy.ndarray of shape (:,3,3) or (3,3)
F : numpy.ndarray of shape (...,3,3)
Deformation gradient.
P : numpy.ndarray of shape (:,3,3) or (3,3)
First Piola-Kirchhoff stress.
Returns
-------
B : numpy.ndarray of shape (...,3,3)
Left Cauchy-Green deformation tensor.
"""
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)
return symmetric(sigma)
return _np.matmul(F,_tensor.transpose(F))
def deviatoric_part(T):
def deformation_Cauchy_Green_right(F):
"""
Return deviatoric part of a tensor.
Calculate right Cauchy-Green deformation tensor.
Parameters
----------
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the deviatoric part is computed.
F : numpy.ndarray of shape (...,3,3)
Deformation gradient.
Returns
-------
C : numpy.ndarray of shape (...,3,3)
Right Cauchy-Green deformation tensor.
"""
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 _np.matmul(_tensor.transpose(F),F)
def eigenvalues(T_sym):
def equivalent_strain_Mises(epsilon):
"""
Return the eigenvalues, i.e. principal components, of a symmetric tensor.
The eigenvalues are sorted in ascending order, each repeated according to
its multiplicity.
Calculate the Mises equivalent of a strain tensor.
Parameters
----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvalues are computed.
epsilon : numpy.ndarray of shape (...,3,3)
Symmetric strain tensor of which the von Mises equivalent is computed.
Returns
-------
epsilon_vM : numpy.ndarray of shape (...)
Von Mises equivalent strain of epsilon.
"""
return _np.linalg.eigvalsh(symmetric(T_sym))
return _equivalent_Mises(epsilon,2.0/3.0)
def eigenvectors(T_sym,RHS=False):
def equivalent_stress_Mises(sigma):
"""
Return eigenvectors of a symmetric tensor.
The eigenvalues are sorted in ascending order of their associated eigenvalues.
Calculate the Mises equivalent of a stress tensor.
Parameters
----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric tensor of which the eigenvectors are computed.
RHS: bool, optional
Enforce right-handed coordinate system. Default is False.
sigma : numpy.ndarray of shape (...,3,3)
Symmetric stress tensor of which the von Mises equivalent is computed.
Returns
-------
sigma_vM : numpy.ndarray of shape (...)
Von Mises equivalent stress of sigma.
"""
(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
else:
v[_np.linalg.det(v) < 0.0,:,2] *= -1.0
return v
def left_stretch(T):
"""
Return the left stretch of a tensor.
Parameters
----------
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the left stretch is computed.
"""
return _polar_decomposition(T,'V')[0]
return _equivalent_Mises(sigma,3.0/2.0)
def maximum_shear(T_sym):
"""
Return the maximum shear component of a symmetric tensor.
Calculate the maximum shear component of a symmetric tensor.
Parameters
----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the maximum shear is computed.
Returns
-------
gamma_max : numpy.ndarray of shape (...)
Maximum shear of T_sym.
"""
w = eigenvalues(T_sym)
return (w[0] - w[2])*0.5 if _np.shape(T_sym) == (3,3) else \
(w[:,0] - w[:,2])*0.5
w = _tensor.eigenvalues(T_sym)
return (w[...,0] - w[...,2])*0.5
def Mises_strain(epsilon):
def rotation(T):
"""
Return the Mises equivalent of a strain tensor.
Calculate the rotational part of a tensor.
Parameters
----------
epsilon : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric strain tensor of which the von Mises equivalent is computed.
"""
return _Mises(epsilon,2.0/3.0)
def Mises_stress(sigma):
"""
Return the Mises equivalent of a stress tensor.
Parameters
----------
sigma : numpy.ndarray of shape (:,3,3) or (3,3)
Symmetric stress tensor of which the von Mises equivalent is computed.
"""
return _Mises(sigma,3.0/2.0)
def PK2(P,F):
"""
Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient.
Parameters
----------
P : numpy.ndarray of shape (...,3,3) or (3,3)
First Piola-Kirchhoff stress.
F : numpy.ndarray of shape (...,3,3) or (3,3)
Deformation gradient.
"""
if _np.shape(F) == _np.shape(P) == (3,3):
S = _np.dot(_np.linalg.inv(F),P)
else:
S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
return symmetric(S)
def right_stretch(T):
"""
Return the right stretch of a tensor.
Parameters
----------
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the right stretch is computed.
"""
return _polar_decomposition(T,'U')[0]
def rotational_part(T):
"""
Return the rotational part of a tensor.
Parameters
----------
T : numpy.ndarray of shape (:,3,3) or (3,3)
T : numpy.ndarray of shape (...,3,3)
Tensor of which the rotational part is computed.
"""
return _polar_decomposition(T,'R')[0]
def spherical_part(T,tensor=False):
"""
Return spherical (hydrostatic) part of a tensor.
Parameters
----------
T : numpy.ndarray of shape (:,3,3) or (3,3)
Tensor of which the hydrostatic part is computed.
tensor : bool, optional
Map spherical part onto identity tensor. Default is false
Returns
-------
R : damask.Rotation of shape (...)
Rotational part of the vector.
"""
if T.shape == (3,3):
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
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 _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0])
def strain_tensor(F,t,m):
def strain(F,t,m):
"""
Return strain tensor calculated from deformation gradient.
Calculate strain tensor (Seth–Hill family).
For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
https://de.wikipedia.org/wiki/Verzerrungstensor
Parameters
----------
F : numpy.ndarray of shape (:,3,3) or (3,3)
F : numpy.ndarray of shape (...,3,3)
Deformation gradient.
t : {‘V’, ‘U’}
Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
Type of the polar decomposition, ‘V’ for left stretch tensor
and ‘U’ for right stretch tensor.
m : float
Order of the strain.
Returns
-------
epsilon : numpy.ndarray of shape (...,3,3)
Strain of F.
"""
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)
w,n = _np.linalg.eigh(deformation_Cauchy_Green_left(F))
elif t == 'U':
C = _np.matmul(transpose(F_),F_)
w,n = _np.linalg.eigh(C)
w,n = _np.linalg.eigh(deformation_Cauchy_Green_right(F))
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.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(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.einsum('...j,...kj,...lj',w**m,n,n) + _np.eye(3))
else:
eps = _np.matmul(n,_np.einsum('ij,ikj->ijk',0.5*_np.log(w),n))
eps = _np.einsum('...j,...kj,...lj',0.5*_np.log(w),n,n)
return eps.reshape(3,3) if _np.shape(F) == (3,3) else \
eps
return eps
def symmetric(T):
def stress_Cauchy(P,F):
"""
Return the symmetrized tensor.
Calculate the Cauchy stress (true stress).
Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
Parameters
----------
T : numpy.ndarray of shape (...,3,3) or (3,3)
Tensor of which the symmetrized values are computed.
P : numpy.ndarray of shape (...,3,3)
First Piola-Kirchhoff stress.
F : numpy.ndarray of shape (...,3,3)
Deformation gradient.
Returns
-------
sigma : numpy.ndarray of shape (...,3,3)
Cauchy stress.
"""
return (T+transpose(T))*0.5
return _tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F))
def transpose(T):
def stress_second_Piola_Kirchhoff(P,F):
"""
Return the transpose of a tensor.
Calculate the second Piola-Kirchhoff stress.
Resulting tensor is symmetrized as the second Piola-Kirchhoff stress
needs to be symmetric.
Parameters
----------
T : numpy.ndarray of shape (...,3,3) or (3,3)
Tensor of which the transpose is computed.
P : numpy.ndarray of shape (...,3,3)
First Piola-Kirchhoff stress.
F : numpy.ndarray of shape (...,3,3)
Deformation gradient.
Returns
-------
S : numpy.ndarray of shape (...,3,3)
Second Piola-Kirchhoff stress.
"""
return T.T if _np.shape(T) == (3,3) else \
_np.swapaxes(T,axis2=-2,axis1=-1)
return _tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P))
def stretch_left(T):
"""
Calculate left stretch of a tensor.
Parameters
----------
T : numpy.ndarray of shape (...,3,3)
Tensor of which the left stretch is computed.
Returns
-------
V : numpy.ndarray of shape (...,3,3)
Left stretch tensor from Polar decomposition of T.
"""
return _polar_decomposition(T,'V')[0]
def stretch_right(T):
"""
Calculate right stretch of a tensor.
Parameters
----------
T : numpy.ndarray of shape (...,3,3)
Tensor of which the right stretch is computed.
Returns
-------
U : numpy.ndarray of shape (...,3,3)
Left stretch tensor from Polar decomposition of T.
"""
return _polar_decomposition(T,'U')[0]
def _polar_decomposition(T,requested):
"""
Singular value decomposition.
Perform singular value decomposition.
Parameters
----------
T : numpy.ndarray of shape (:,3,3) or (3,3)
T : numpy.ndarray of shape (...,3,3)
Tensor of which the singular values are computed.
requested : iterable of str
Requested outputs: ‘R’ for the rotation tensor,
‘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, _, vh = _np.linalg.svd(T)
R = _np.einsum('...ij,...jk',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.einsum('...ij,...kj',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.einsum('...ji,...jk',R,T))
if len(output) == 0:
raise ValueError('output needs to be out of V, R, U')
return tuple(output)
def _Mises(T_sym,s):
def _equivalent_Mises(T_sym,s):
"""
Base equation for Mises equivalent of a stres or strain tensor.
Base equation for Mises equivalent of a stress or strain tensor.
Parameters
----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3)
T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the von Mises equivalent is computed.
s : float
Scaling factor (2/3 for strain, 3/2 for stress).
"""
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))
d = _tensor.deviatoric(T_sym)
return _np.sqrt(s*_np.sum(d**2.0,axis=(-1,-2)))

2
python/damask/seeds.py
View File

@ -1,3 +1,5 @@
"""Functionality for generation of seed points for Voronoi or Laguerre tessellation."""
from scipy import spatial as _spatial
import numpy as _np

131
python/damask/tensor.py Normal file
View File

@ -0,0 +1,131 @@
"""
Tensor operations.
Notes
-----
This is not a tensor class, but a collection of routines
to operate on numpy.ndarrays of shape (...,3,3).
"""
import numpy as _np
def deviatoric(T):
"""
Calculate deviatoric part of a tensor.
Parameters
----------
T : numpy.ndarray of shape (...,3,3)
Tensor of which the deviatoric part is computed.
Returns
-------
T' : numpy.ndarray of shape (...,3,3)
Deviatoric part of T.
"""
return T - spherical(T,tensor=True)
def eigenvalues(T_sym):
"""
Eigenvalues, i.e. principal components, of a symmetric tensor.
Parameters
----------
T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the eigenvalues are computed.
Returns
-------
lambda : numpy.ndarray of shape (...,3)
Eigenvalues of T_sym sorted in ascending order, each repeated
according to its multiplicity.
"""
return _np.linalg.eigvalsh(symmetric(T_sym))
def eigenvectors(T_sym,RHS=False):
"""
Eigenvectors of a symmetric tensor.
Parameters
----------
T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the eigenvectors are computed.
RHS: bool, optional
Enforce right-handed coordinate system. Defaults to False.
Returns
-------
x : numpy.ndarray of shape (...,3,3)
Eigenvectors of T_sym sorted in ascending order of their
associated eigenvalues.
"""
(u,v) = _np.linalg.eigh(symmetric(T_sym))
if RHS:
v[_np.linalg.det(v) < 0.0,:,2] *= -1.0
return v
def spherical(T,tensor=True):
"""
Calculate spherical part of a tensor.
Parameters
----------
T : numpy.ndarray of shape (...,3,3)
Tensor of which the spherical part is computed.
tensor : bool, optional
Map spherical part onto identity tensor. Defaults to True.
Returns
-------
p : numpy.ndarray of shape (...,3,3)
unless tensor == False: shape (...,)
Spherical part of tensor T. p is an isotropic tensor.
"""
sph = _np.trace(T,axis2=-2,axis1=-1)/3.0
return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph
def symmetric(T):
"""
Symmetrize tensor.
Parameters
----------
T : numpy.ndarray of shape (...,3,3)
Tensor of which the symmetrized values are computed.
Returns
-------
T_sym : numpy.ndarray of shape (...,3,3)
Symmetrized tensor T.
"""
return (T+transpose(T))*0.5
def transpose(T):
"""
Transpose tensor.
Parameters
----------
T : numpy.ndarray of shape (...,3,3)
Tensor of which the transpose is computed.
Returns
-------
T.T : numpy.ndarray of shape (...,3,3)
Transpose of tensor T.
"""
return _np.swapaxes(T,axis2=-2,axis1=-1)

2
python/tests/conftest.py
View File

@ -82,7 +82,7 @@ def set_of_quaternions():
return qu
n = 1100
n = 600
scatter=1.e-2
specials = np.array([
[1.0, 0.0, 0.0, 0.0],

4
python/tests/test_Geom.py
View File

@ -55,7 +55,7 @@ class TestGeom:
def test_invalid_vtr(self,tmp_path):
v = VTK.from_rectilinear_grid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0)
v.save(tmp_path/'no_materialpoint.vtr')
v.save(tmp_path/'no_materialpoint.vtr',parallel=False)
with pytest.raises(ValueError):
Geom.load(tmp_path/'no_materialpoint.vtr')
@ -229,7 +229,7 @@ class TestGeom:
@pytest.mark.parametrize('Eulers',[[32.0,68.0,21.0],
[0.0,32.0,240.0]])
def test_rotate(self,default,update,reference_dir,Eulers):
modified = default.rotate(Rotation.from_Eulers(Eulers,degrees=True))
modified = default.rotate(Rotation.from_Euler_angles(Eulers,degrees=True))
tag = f'Eulers_{util.srepr(Eulers,"-")}'
reference = reference_dir/f'rotate_{tag}.vtr'
if update: modified.save(reference)

22
python/tests/test_Orientation.py
View File

@ -16,7 +16,7 @@ def reference_dir(reference_dir_base):
@pytest.fixture
def set_of_rodrigues(set_of_quaternions):
return Rotation(set_of_quaternions).as_Rodrigues()[:200]
return Rotation(set_of_quaternions).as_Rodrigues_vector()
class TestOrientation:
@ -90,8 +90,8 @@ class TestOrientation:
assert np.all(Orientation.from_quaternion(q=np.array([1,0,0,0]),lattice='triclinic').as_matrix()
== np.eye(3))
def test_from_Eulers(self):
assert np.all(Orientation.from_Eulers(phi=np.zeros(3),lattice='triclinic').as_matrix()
def test_from_Euler_angles(self):
assert np.all(Orientation.from_Euler_angles(phi=np.zeros(3),lattice='triclinic').as_matrix()
== np.eye(3))
def test_from_axis_angle(self):
@ -106,8 +106,8 @@ class TestOrientation:
assert np.all(Orientation.from_matrix(R=np.eye(3),lattice='triclinic').as_matrix()
== np.eye(3))
def test_from_Rodrigues(self):
assert np.all(Orientation.from_Rodrigues(rho=np.array([0,0,1,0]),lattice='triclinic').as_matrix()
def test_from_Rodrigues_vector(self):
assert np.all(Orientation.from_Rodrigues_vector(rho=np.array([0,0,1,0]),lattice='triclinic').as_matrix()
== np.eye(3))
def test_from_homochoric(self):
@ -208,7 +208,7 @@ class TestOrientation:
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
def test_disorientation360(self,lattice):
o_1 = Orientation(Rotation(),lattice)
o_2 = Orientation.from_Eulers(lattice=lattice,phi=[360,0,0],degrees=True)
o_2 = Orientation.from_Euler_angles(lattice=lattice,phi=[360,0,0],degrees=True)
assert np.allclose((o_1.disorientation(o_2)).as_matrix(),np.eye(3))
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@ -275,16 +275,16 @@ class TestOrientation:
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
def test_in_FZ_vectorization(self,set_of_rodrigues,lattice):
result = Orientation.from_Rodrigues(rho=set_of_rodrigues.reshape((50,4,-1)),lattice=lattice).in_FZ.reshape(-1)
result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),lattice=lattice).in_FZ.reshape(-1)
for r,rho in zip(result,set_of_rodrigues[:len(result)]):
assert r == Orientation.from_Rodrigues(rho=rho,lattice=lattice).in_FZ
assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_FZ
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
def test_in_disorientation_FZ_vectorization(self,set_of_rodrigues,lattice):
result = Orientation.from_Rodrigues(rho=set_of_rodrigues.reshape((50,4,-1)),
result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),
lattice=lattice).in_disorientation_FZ.reshape(-1)
for r,rho in zip(result,set_of_rodrigues[:len(result)]):
assert r == Orientation.from_Rodrigues(rho=rho,lattice=lattice).in_disorientation_FZ
assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_disorientation_FZ
@pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@ -417,7 +417,7 @@ class TestOrientation:
def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = reference_dir/f'{lattice}_{model}.txt'
o = Orientation(lattice=lattice)
eu = o.related(model).as_Eulers(degrees=True)
eu = o.related(model).as_Euler_angles(degrees=True)
if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)])
Table(eu,{'Eulers':(3,)})\

69
python/tests/test_Result.py
View File

@ -11,6 +11,7 @@ import h5py
from damask import Result
from damask import Rotation
from damask import Orientation
from damask import tensor
from damask import mechanics
from damask import grid_filters
@ -120,13 +121,13 @@ class TestResult:
in_file = default.read_dataset(loc['x'],0)
assert np.allclose(in_memory,in_file)
def test_add_Cauchy(self,default):
default.add_Cauchy('P','F')
def test_add_stress_Cauchy(self,default):
default.add_stress_Cauchy('P','F')
loc = {'F': default.get_dataset_location('F'),
'P': default.get_dataset_location('P'),
'sigma':default.get_dataset_location('sigma')}
in_memory = mechanics.Cauchy(default.read_dataset(loc['P'],0),
default.read_dataset(loc['F'],0))
in_memory = mechanics.stress_Cauchy(default.read_dataset(loc['P'],0),
default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['sigma'],0)
assert np.allclose(in_memory,in_file)
@ -142,27 +143,27 @@ class TestResult:
default.add_deviator('P')
loc = {'P' :default.get_dataset_location('P'),
's_P':default.get_dataset_location('s_P')}
in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0))
in_memory = tensor.deviatoric(default.read_dataset(loc['P'],0))
in_file = default.read_dataset(loc['s_P'],0)
assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('eigenvalue,function',[('max',np.amax),('min',np.amin)])
def test_add_eigenvalue(self,default,eigenvalue,function):
default.add_Cauchy('P','F')
default.add_stress_Cauchy('P','F')
default.add_eigenvalue('sigma',eigenvalue)
loc = {'sigma' :default.get_dataset_location('sigma'),
'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')}
in_memory = function(mechanics.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True)
in_memory = function(tensor.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True)
in_file = default.read_dataset(loc['lambda'],0)
assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('eigenvalue,idx',[('max',2),('mid',1),('min',0)])
def test_add_eigenvector(self,default,eigenvalue,idx):
default.add_Cauchy('P','F')
default.add_stress_Cauchy('P','F')
default.add_eigenvector('sigma',eigenvalue)
loc = {'sigma' :default.get_dataset_location('sigma'),
'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')}
in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx]
in_memory = tensor.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx]
in_file = default.read_dataset(loc['v(sigma)'],0)
assert np.allclose(in_memory,in_file)
@ -182,7 +183,7 @@ class TestResult:
assert np.allclose(in_memory,in_file)
def test_add_maximum_shear(self,default):
default.add_Cauchy('P','F')
default.add_stress_Cauchy('P','F')
default.add_maximum_shear('sigma')
loc = {'sigma' :default.get_dataset_location('sigma'),
'max_shear(sigma)':default.get_dataset_location('max_shear(sigma)')}
@ -193,36 +194,36 @@ class TestResult:
def test_add_Mises_strain(self,default):
t = ['V','U'][np.random.randint(0,2)]
m = np.random.random()*2.0 - 1.0
default.add_strain_tensor('F',t,m)
default.add_strain('F',t,m)
label = f'epsilon_{t}^{m}(F)'
default.add_Mises(label)
default.add_equivalent_Mises(label)
loc = {label :default.get_dataset_location(label),
label+'_vM':default.get_dataset_location(label+'_vM')}
in_memory = mechanics.Mises_strain(default.read_dataset(loc[label],0)).reshape(-1,1)
in_memory = mechanics.equivalent_strain_Mises(default.read_dataset(loc[label],0)).reshape(-1,1)
in_file = default.read_dataset(loc[label+'_vM'],0)
assert np.allclose(in_memory,in_file)
def test_add_Mises_stress(self,default):
default.add_Cauchy('P','F')
default.add_Mises('sigma')
default.add_stress_Cauchy('P','F')
default.add_equivalent_Mises('sigma')
loc = {'sigma' :default.get_dataset_location('sigma'),
'sigma_vM':default.get_dataset_location('sigma_vM')}
in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
in_memory = mechanics.equivalent_stress_Mises(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
in_file = default.read_dataset(loc['sigma_vM'],0)
assert np.allclose(in_memory,in_file)
def test_add_Mises_invalid(self,default):
default.add_Cauchy('P','F')
default.add_stress_Cauchy('P','F')
default.add_calculation('sigma_y','#sigma#',unit='y')
default.add_Mises('sigma_y')
default.add_equivalent_Mises('sigma_y')
assert default.get_dataset_location('sigma_y_vM') == []
def test_add_Mises_stress_strain(self,default):
default.add_Cauchy('P','F')
default.add_stress_Cauchy('P','F')
default.add_calculation('sigma_y','#sigma#',unit='y')
default.add_calculation('sigma_x','#sigma#',unit='x')
default.add_Mises('sigma_y',kind='strain')
default.add_Mises('sigma_x',kind='stress')
default.add_equivalent_Mises('sigma_y',kind='strain')
default.add_equivalent_Mises('sigma_x',kind='stress')
loc = {'y' :default.get_dataset_location('sigma_y_vM'),
'x' :default.get_dataset_location('sigma_x_vM')}
assert not np.allclose(default.read_dataset(loc['y'],0),default.read_dataset(loc['x'],0))
@ -235,13 +236,13 @@ class TestResult:
in_file = default.read_dataset(loc['|F|_1'],0)
assert np.allclose(in_memory,in_file)
def test_add_PK2(self,default):
default.add_PK2('P','F')
def test_add_stress_second_Piola_Kirchhoff(self,default):
default.add_stress_second_Piola_Kirchhoff('P','F')
loc = {'F':default.get_dataset_location('F'),
'P':default.get_dataset_location('P'),
'S':default.get_dataset_location('S')}
in_memory = mechanics.PK2(default.read_dataset(loc['P'],0),
default.read_dataset(loc['F'],0))
in_memory = mechanics.stress_second_Piola_Kirchhoff(default.read_dataset(loc['P'],0),
default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['S'],0)
assert np.allclose(in_memory,in_file)
@ -260,11 +261,11 @@ class TestResult:
in_file = default.read_dataset(loc['pole'])
assert np.allclose(in_memory,in_file)
def test_add_rotational_part(self,default):
default.add_rotational_part('F')
def test_add_rotation(self,default):
default.add_rotation('F')
loc = {'F': default.get_dataset_location('F'),
'R(F)': default.get_dataset_location('R(F)')}
in_memory = mechanics.rotational_part(default.read_dataset(loc['F'],0))
in_memory = mechanics.rotation(default.read_dataset(loc['F'],0)).as_matrix()
in_file = default.read_dataset(loc['R(F)'],0)
assert np.allclose(in_memory,in_file)
@ -272,18 +273,18 @@ class TestResult:
default.add_spherical('P')
loc = {'P': default.get_dataset_location('P'),
'p_P': default.get_dataset_location('p_P')}
in_memory = mechanics.spherical_part(default.read_dataset(loc['P'],0)).reshape(-1,1)
in_memory = tensor.spherical(default.read_dataset(loc['P'],0),False).reshape(-1,1)
in_file = default.read_dataset(loc['p_P'],0)
assert np.allclose(in_memory,in_file)
def test_add_strain(self,default):
t = ['V','U'][np.random.randint(0,2)]
m = np.random.random()*2.0 - 1.0
default.add_strain_tensor('F',t,m)
default.add_strain('F',t,m)
label = f'epsilon_{t}^{m}(F)'
loc = {'F': default.get_dataset_location('F'),
label: default.get_dataset_location(label)}
in_memory = mechanics.strain_tensor(default.read_dataset(loc['F'],0),t,m)
in_memory = mechanics.strain(default.read_dataset(loc['F'],0),t,m)
in_file = default.read_dataset(loc[label],0)
assert np.allclose(in_memory,in_file)
@ -291,7 +292,7 @@ class TestResult:
default.add_stretch_tensor('F','U')
loc = {'F': default.get_dataset_location('F'),
'U(F)': default.get_dataset_location('U(F)')}
in_memory = mechanics.right_stretch(default.read_dataset(loc['F'],0))
in_memory = mechanics.stretch_right(default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['U(F)'],0)
assert np.allclose(in_memory,in_file)
@ -299,7 +300,7 @@ class TestResult:
default.add_stretch_tensor('F','V')
loc = {'F': default.get_dataset_location('F'),
'V(F)': default.get_dataset_location('V(F)')}
in_memory = mechanics.left_stretch(default.read_dataset(loc['F'],0))
in_memory = mechanics.stretch_left(default.read_dataset(loc['F'],0))
in_file = default.read_dataset(loc['V(F)'],0)
assert np.allclose(in_memory,in_file)
@ -311,7 +312,7 @@ class TestResult:
def test_add_overwrite(self,default,overwrite):
default.pick('times',default.times_in_range(0,np.inf)[-1])
default.add_Cauchy()
default.add_stress_Cauchy()
loc = default.get_dataset_location('sigma')
with h5py.File(default.fname,'r') as f:
# h5py3 compatibility

66
python/tests/test_Rotation.py
View File

@ -522,7 +522,7 @@ class TestRotation:
def test_Eulers_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from Euler angles and back."""
for rot in set_of_rotations:
m = rot.as_Eulers()
m = rot.as_Euler_angles()
o = backward(forward(m))
u = np.array([np.pi*2,np.pi,np.pi*2])
ok = np.allclose(m,o,atol=atol)
@ -559,7 +559,7 @@ class TestRotation:
"""Ensure invariance of conversion from Rodrigues-Frank vector and back."""
cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in set_of_rotations:
m = rot.as_Rodrigues()
m = rot.as_Rodrigues_vector()
o = backward(forward(m))
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol)
@ -626,7 +626,7 @@ class TestRotation:
(Rotation._eu2ro,eu2ro)])
def test_Eulers_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for Euler angles against single point calculation."""
eu = np.array([rot.as_Eulers() for rot in set_of_rotations])
eu = np.array([rot.as_Euler_angles() for rot in set_of_rotations])
vectorized(eu.reshape(eu.shape[0]//2,-1,3))
co = vectorized(eu)
for e,c in zip(eu,co):
@ -649,7 +649,7 @@ class TestRotation:
(Rotation._ro2ho,ro2ho)])
def test_Rodrigues_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for Rodrigues-Frank vector against single point calculation."""
ro = np.array([rot.as_Rodrigues() for rot in set_of_rotations])
ro = np.array([rot.as_Rodrigues_vector() for rot in set_of_rotations])
vectorized(ro.reshape(ro.shape[0]//2,-1,4))
co = vectorized(ro)
for r,c in zip(ro,co):
@ -682,12 +682,15 @@ class TestRotation:
for v,o in zip(vec,ori):
assert np.allclose(func(v,normalize=True).as_quaternion(),o.as_quaternion())
def test_invalid_init(self):
with pytest.raises(TypeError):
Rotation(np.ones(3))
@pytest.mark.parametrize('degrees',[True,False])
def test_Eulers(self,set_of_rotations,degrees):
for rot in set_of_rotations:
m = rot.as_quaternion()
o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion()
o = Rotation.from_Euler_angles(rot.as_Euler_angles(degrees),degrees).as_quaternion()
ok = np.allclose(m,o,atol=atol)
if np.isclose(rot.as_quaternion()[0],0.0,atol=atol):
ok |= np.allclose(m*-1.,o,atol=atol)
@ -699,8 +702,8 @@ class TestRotation:
def test_axis_angle(self,set_of_rotations,degrees,normalize,P):
c = np.array([P*-1,P*-1,P*-1,1.])
for rot in set_of_rotations:
m = rot.as_Eulers()
o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalize,P).as_Eulers()
m = rot.as_Euler_angles()
o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalize,P).as_Euler_angles()
u = np.array([np.pi*2,np.pi,np.pi*2])
ok = np.allclose(m,o,atol=atol)
ok |= np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
@ -726,7 +729,7 @@ class TestRotation:
c = np.array([P*-1,P*-1,P*-1,1.])
for rot in set_of_rotations:
m = rot.as_matrix()
o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalize,P).as_matrix()
o = Rotation.from_Rodrigues_vector(rot.as_Rodrigues_vector()*c,normalize,P).as_matrix()
ok = np.allclose(m,o,atol=atol)
assert ok and np.isclose(np.linalg.det(o),1.0), f'{m},{o}'
@ -734,8 +737,8 @@ class TestRotation:
def test_homochoric(self,set_of_rotations,P):
cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in set_of_rotations:
m = rot.as_Rodrigues()
o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues()
m = rot.as_Rodrigues_vector()
o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues_vector()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol)
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}'
@ -820,10 +823,10 @@ class TestRotation:
assert r == r.flatten(order).reshape(shape,order)
@pytest.mark.parametrize('function',[Rotation.from_quaternion,
Rotation.from_Eulers,
Rotation.from_Euler_angles,
Rotation.from_axis_angle,
Rotation.from_matrix,
Rotation.from_Rodrigues,
Rotation.from_Rodrigues_vector,
Rotation.from_homochoric,
Rotation.from_cubochoric])
def test_invalid_shape(self,function):
@ -831,9 +834,16 @@ class TestRotation:
with pytest.raises(ValueError):
function(invalid_shape)
def test_invalid_shape_parallel(self):
invalid_a = np.random.random(np.random.randint(8,32,(3)))
invalid_b = np.random.random(np.random.randint(8,32,(3)))
with pytest.raises(ValueError):
Rotation.from_parallel(invalid_a,invalid_b)
@pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'),
(Rotation.from_axis_angle,'as_axis_angle'),
(Rotation.from_Rodrigues, 'as_Rodrigues'),
(Rotation.from_Rodrigues_vector, 'as_Rodrigues_vector'),
(Rotation.from_homochoric,'as_homochoric'),
(Rotation.from_cubochoric,'as_cubochoric')])
def test_invalid_P(self,fr,to):
@ -852,16 +862,16 @@ class TestRotation:
@pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
(Rotation.from_quaternion, np.array([1,1,1,0])),
(Rotation.from_Eulers, np.array([1,4,0])),
(Rotation.from_axis_angle, np.array([1,0,0,4])),
(Rotation.from_axis_angle, np.array([1,1,0,1])),
(Rotation.from_matrix, np.random.rand(3,3)),
(Rotation.from_matrix, np.array([[1,1,0],[1,2,0],[0,0,1]])),
(Rotation.from_Rodrigues, np.array([1,0,0,-1])),
(Rotation.from_Rodrigues, np.array([1,1,0,1])),
(Rotation.from_homochoric, np.array([2,2,2])),
(Rotation.from_cubochoric, np.array([1.1,0,0])) ])
(Rotation.from_quaternion, np.array([1,1,1,0])),
(Rotation.from_Euler_angles, np.array([1,4,0])),
(Rotation.from_axis_angle, np.array([1,0,0,4])),
(Rotation.from_axis_angle, np.array([1,1,0,1])),
(Rotation.from_matrix, np.random.rand(3,3)),
(Rotation.from_matrix, np.array([[1,1,0],[1,2,0],[0,0,1]])),
(Rotation.from_Rodrigues_vector, np.array([1,0,0,-1])),
(Rotation.from_Rodrigues_vector, np.array([1,1,0,1])),
(Rotation.from_homochoric, np.array([2,2,2])),
(Rotation.from_cubochoric, np.array([1.1,0,0])) ])
def test_invalid_value(self,function,invalid):
with pytest.raises(ValueError):
function(invalid)
@ -905,8 +915,8 @@ class TestRotation:
def test_rotate_360deg(self,data):
phi_1 = np.random.random() * np.pi
phi_2 = 2*np.pi - phi_1
R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.]))
R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2]))
R_1 = Rotation.from_Euler_angles(np.array([phi_1,0.,0.]))
R_2 = Rotation.from_Euler_angles(np.array([0.,0.,phi_2]))
assert np.allclose(data,R_2@(R_1@data))
@pytest.mark.parametrize('pwr',[-10,0,1,2.5,np.pi,np.random.random()])
@ -952,7 +962,7 @@ class TestRotation:
def test_misorientation360(self):
R_1 = Rotation()
R_2 = Rotation.from_Eulers([360,0,0],degrees=True)
R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True)
assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3))
@pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120])
@ -1015,7 +1025,7 @@ class TestRotation:
Eulers = grid_filters.cell_coord0(steps,limits)
Eulers = np.radians(Eulers) if not degrees else Eulers
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Eulers(True)
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
if fractions: assert np.sqrt(((weights_r - weights) ** 2).mean()) < 4
@ -1033,7 +1043,7 @@ class TestRotation:
Eulers = grid_filters.node_coord0(steps,limits)[:-1,:-1,:-1]
Eulers = np.radians(Eulers) if not degrees else Eulers
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Eulers(True)
Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Euler_angles(True)
weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights)
assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5

15
python/tests/test_VTK.py
View File

@ -100,12 +100,17 @@ class TestVTK:
v = VTK.from_poly_data(points)
v.save(tmp_path/fname)
@pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk', None),
('this_file_does_not_exist.vtk','vtk'),
('this_file_does_not_exist.vtx', None)])
def test_invalid_dataset_type(self,name,dataset_type):
@pytest.mark.parametrize('fname,dataset_type',[('a_file.vtk', None),
('a_file.vtk','vtk'),
('a_file.vtx', None)])
def test_invalid_dataset_type(self,tmp_path,fname,dataset_type):
open(tmp_path/fname,'a').close()
with pytest.raises(TypeError):
VTK.load(name,dataset_type)
VTK.load(tmp_path/fname,dataset_type)
def test_file_not_found(self):
with pytest.raises(FileNotFoundError):
VTK.load('/dev/null')
def test_add_extension(self,tmp_path,default):
default.save(tmp_path/'default.txt',parallel=False)

260
python/tests/test_mechanics.py
View File

@ -1,88 +1,177 @@
import pytest
import numpy as np
from damask import tensor
from damask import mechanics
from damask import Rotation
def stress_Cauchy(P,F):
sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T)
return symmetric(sigma)
def eigenvalues(T_sym):
return np.linalg.eigvalsh(symmetric(T_sym))
def maximum_shear(T_sym):
w = eigenvalues(T_sym)
return (w[0] - w[2])*0.5
def equivalent_strain_Mises(epsilon):
return equivalent_Mises(epsilon,2.0/3.0)
def equivalent_stress_Mises(sigma):
return equivalent_Mises(sigma,3.0/2.0)
def stress_second_Piola_Kirchhoff(P,F):
S = np.dot(np.linalg.inv(F),P)
return symmetric(S)
def rotation(T):
return polar_decomposition(T,'R')[0]
def strain(F,t,m):
if t == 'V':
B = np.matmul(F,F.T)
w,n = np.linalg.eigh(B)
elif t == 'U':
C = np.matmul(F.T,F)
w,n = np.linalg.eigh(C)
if m > 0.0:
eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('j,kj->jk',w**m,n))
- np.eye(3))
elif m < 0.0:
eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('j,kj->jk',w**m,n))
+ np.eye(3))
else:
eps = np.matmul(n,np.einsum('j,kj->jk',0.5*np.log(w),n))
return eps
def stretch_left(T):
return polar_decomposition(T,'V')[0]
def stretch_right(T):
return polar_decomposition(T,'U')[0]
def symmetric(T):
return (T+T.T)*0.5
def polar_decomposition(T,requested):
u, s, vh = np.linalg.svd(T)
R = np.dot(u,vh)
output = []
if 'R' in requested:
output.append(R)
if 'V' in requested:
output.append(np.dot(T,R.T))
if 'U' in requested:
output.append(np.dot(R.T,T))
return tuple(output)
def equivalent_Mises(T_sym,s):
return np.sqrt(s*(np.sum(deviatoric(T_sym)**2.0)))
def deviatoric(T):
return T - np.eye(3)*np.trace(T)/3.0
class TestMechanics:
n = 1000
c = np.random.randint(n)
@pytest.mark.parametrize('vectorized,single',[(mechanics.maximum_shear, maximum_shear),
(mechanics.equivalent_stress_Mises, equivalent_stress_Mises),
(mechanics.equivalent_strain_Mises, equivalent_strain_Mises),
(mechanics.stretch_left, stretch_left),
(mechanics.stretch_right, stretch_right),
])
def test_vectorize_1_arg(self,vectorized,single):
epsilon = np.random.rand(self.n,3,3)
epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3))
for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)):
assert np.allclose(single(epsilon[i]),v)
@pytest.mark.parametrize('function',[mechanics.deviatoric_part,
mechanics.eigenvalues,
mechanics.eigenvectors,
mechanics.left_stretch,
mechanics.maximum_shear,
mechanics.Mises_strain,
mechanics.Mises_stress,
mechanics.right_stretch,
mechanics.rotational_part,
mechanics.spherical_part,
mechanics.symmetric,
mechanics.transpose,
])
def test_vectorize_1_arg(self,function):
epsilon = np.random.rand(self.n,3,3)
assert np.allclose(function(epsilon)[self.c],function(epsilon[self.c]))
def test_vectorize_rotation(self):
epsilon = Rotation.from_random(self.n).as_matrix()
epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3))
for i,v in enumerate(np.reshape(mechanics.rotation(epsilon_vec).as_matrix(),
mechanics.rotation(epsilon).as_matrix().shape)):
assert np.allclose(rotation(epsilon[i]),v)
@pytest.mark.parametrize('function',[mechanics.Cauchy,
mechanics.PK2,
])
def test_vectorize_2_arg(self,function):
P = np.random.rand(self.n,3,3)
F = np.random.rand(self.n,3,3)
assert np.allclose(function(P,F)[self.c],function(P[self.c],F[self.c]))
def test_vectorize_strain_tensor(self):
F = np.random.rand(self.n,3,3)
t = ['V','U'][np.random.randint(0,2)]
m = np.random.random()*10. -5.0
assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c],
mechanics.strain_tensor(F[self.c],t,m))
@pytest.mark.parametrize('vectorized,single',[(mechanics.stress_Cauchy, stress_Cauchy),
(mechanics.stress_second_Piola_Kirchhoff, stress_second_Piola_Kirchhoff)
])
def test_vectorize_2_arg(self,vectorized,single):
P = np.random.rand(self.n,3,3)
F = np.random.rand(self.n,3,3)
P_vec = np.reshape(P,(self.n//10,10,3,3))
F_vec = np.reshape(F,(self.n//10,10,3,3))
for i,v in enumerate(np.reshape(vectorized(P_vec,F_vec),vectorized(P,F).shape)):
assert np.allclose(single(P[i],F[i]),v)
@pytest.mark.parametrize('function',[mechanics.Cauchy,
mechanics.PK2,
@pytest.mark.parametrize('vectorized,single',[(mechanics.strain,strain)])
def test_vectorize_strain(self,vectorized,single):
F = np.random.rand(self.n,3,3)
F_vec = np.reshape(F,(self.n//10,10,3,3))
t = ['V','U'][np.random.randint(0,2)]
m = np.random.random()*10.0 -5.0
for i,v in enumerate(np.reshape(vectorized(F_vec,t,m),vectorized(F,t,m).shape)):
assert np.allclose(single(F[i],t,m),v)
@pytest.mark.parametrize('function',[mechanics.stress_Cauchy,
mechanics.stress_second_Piola_Kirchhoff,
])
def test_stress_measures(self,function):
"""Ensure that all stress measures are equivalent for no deformation."""
P = np.random.rand(self.n,3,3)
assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),mechanics.symmetric(P))
def test_deviatoric_part(self):
I_n = np.broadcast_to(np.eye(3),(self.n,3,3))
r = np.logical_not(I_n)*np.random.rand(self.n,3,3)
assert np.allclose(mechanics.deviatoric_part(I_n+r),r)
assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),tensor.symmetric(P))
def test_polar_decomposition(self):
"""F = RU = VR."""
F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3)
R = mechanics.rotational_part(F)
V = mechanics.left_stretch(F)
U = mechanics.right_stretch(F)
R = mechanics.rotation(F).as_matrix()
V = mechanics.stretch_left(F)
U = mechanics.stretch_right(F)
assert np.allclose(np.matmul(R,U),
np.matmul(V,R))
@pytest.mark.parametrize('m',[0.0,np.random.random()*10.,np.random.random()*-10.])
def test_strain_tensor_no_rotation(self,m):
def test_strain_no_rotation(self,m):
"""Ensure that left and right stretch give same results for no rotation."""
F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3)
assert np.allclose(mechanics.strain_tensor(F,'U',m),
mechanics.strain_tensor(F,'V',m))
assert np.allclose(mechanics.strain(F,'U',m),
mechanics.strain(F,'V',m))
@pytest.mark.parametrize('m',[0.0,np.random.random()*2.5,np.random.random()*-2.5])
def test_strain_tensor_rotation_equivalence(self,m):
def test_strain_rotation_equivalence(self,m):
"""Ensure that left and right strain differ only by a rotation."""
F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.rand(self.n,3,3)*0.5 - 0.25)
assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
np.linalg.det(mechanics.strain_tensor(F,'V',m)))
assert np.allclose(np.linalg.det(mechanics.strain(F,'U',m)),
np.linalg.det(mechanics.strain(F,'V',m)))
@pytest.mark.parametrize('m',[0.0,np.random.random(),np.random.random()*-1.])
@pytest.mark.parametrize('t',['V','U'])
def test_strain_tensor_rotation(self,m,t):
def test_strain_rotation(self,m,t):
"""Ensure that pure rotation results in no strain."""
F = mechanics.rotational_part(np.random.rand(self.n,3,3))
assert np.allclose(mechanics.strain_tensor(F,t,m),
F = Rotation.from_random(self.n).as_matrix()
assert np.allclose(mechanics.strain(F,t,m),
0.0)
def test_rotation_determinant(self):
@ -92,82 +181,37 @@ class TestMechanics:
Should be +1, but random F might contain a reflection.
"""
x = np.random.rand(self.n,3,3)
assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))),
assert np.allclose(np.abs(np.linalg.det(mechanics._polar_decomposition(x,'R')[0])),
1.0)
def test_spherical_deviatoric_part(self):
"""Ensure that full tensor is sum of spherical and deviatoric part."""
x = np.random.rand(self.n,3,3)
sph = mechanics.spherical_part(x,True)
assert np.allclose(sph + mechanics.deviatoric_part(x),
x)
def test_deviatoric_Mises(self):
"""Ensure that Mises equivalent stress depends only on deviatoric part."""
x = np.random.rand(self.n,3,3)
full = mechanics.Mises_stress(x)
dev = mechanics.Mises_stress(mechanics.deviatoric_part(x))
full = mechanics.equivalent_stress_Mises(x)
dev = mechanics.equivalent_stress_Mises(tensor.deviatoric(x))
assert np.allclose(full,
dev)
def test_spherical_mapping(self):
"""Ensure that mapping to tensor is correct."""
@pytest.mark.parametrize('Mises_equivalent',[mechanics.equivalent_strain_Mises,
mechanics.equivalent_stress_Mises])
def test_spherical_Mises(self,Mises_equivalent):
"""Ensure that Mises equivalent strain/stress of spherical strain is 0."""
x = np.random.rand(self.n,3,3)
tensor = mechanics.spherical_part(x,True)
scalar = mechanics.spherical_part(x)
assert np.allclose(np.linalg.det(tensor),
scalar**3.0)
def test_spherical_Mises(self):
"""Ensure that Mises equivalent strrain of spherical strain is 0."""
x = np.random.rand(self.n,3,3)
sph = mechanics.spherical_part(x,True)
assert np.allclose(mechanics.Mises_strain(sph),
assert np.allclose(Mises_equivalent(tensor.spherical(x,True)),
0.0)
def test_symmetric(self):
"""Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
x = np.random.rand(self.n,3,3)
assert np.allclose(mechanics.symmetric(x)*2.0,
mechanics.transpose(x)+x)
def test_transpose(self):
"""Ensure that a symmetric tensor equals its transpose."""
x = mechanics.symmetric(np.random.rand(self.n,3,3))
assert np.allclose(mechanics.transpose(x),
x)
def test_Mises(self):
"""Ensure that equivalent stress is 3/2 of equivalent strain."""
x = np.random.rand(self.n,3,3)
assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
assert np.allclose(mechanics.equivalent_stress_Mises(x)/mechanics.equivalent_strain_Mises(x),
1.5)
def test_eigenvalues(self):
"""Ensure that the characteristic polynomial can be solved."""
A = mechanics.symmetric(np.random.rand(self.n,3,3))
lambd = mechanics.eigenvalues(A)
s = np.random.randint(self.n)
for i in range(3):
assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0)
def test_eigenvalues_and_vectors(self):
"""Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial."""
A = mechanics.symmetric(np.random.rand(self.n,3,3))
lambd = mechanics.eigenvalues(A)
x = mechanics.eigenvectors(A)
s = np.random.randint(self.n)
for i in range(3):
assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0)
def test_eigenvectors_RHS(self):
"""Ensure that RHS coordinate system does only change sign of determinant."""
A = mechanics.symmetric(np.random.rand(self.n,3,3))
LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False))
RHS = np.linalg.det(mechanics.eigenvectors(A,RHS=True))
assert np.allclose(np.abs(LRHS),RHS)
def test_spherical_no_shear(self):
"""Ensure that sherical stress has max shear of 0.0."""
A = mechanics.spherical_part(mechanics.symmetric(np.random.rand(self.n,3,3)),True)
A = tensor.spherical(tensor.symmetric(np.random.rand(self.n,3,3)),True)
assert np.allclose(mechanics.maximum_shear(A),0.0)
def test_invalid_decomposition(self):
with pytest.raises(ValueError):
mechanics._polar_decomposition(np.random.rand(10,3,3),'A')

99
python/tests/test_tensor.py Normal file
View File

@ -0,0 +1,99 @@
import pytest
import numpy as np
from damask import tensor
def deviatoric(T):
return T - spherical(T)
def eigenvalues(T_sym):
return np.linalg.eigvalsh(symmetric(T_sym))
def eigenvectors(T_sym,RHS=False):
(u,v) = np.linalg.eigh(symmetric(T_sym))
if RHS:
if np.linalg.det(v) < 0.0: v[:,2] *= -1.0
return v
def symmetric(T):
return (T+transpose(T))*0.5
def transpose(T):
return T.T
def spherical(T,tensor=True):
sph = np.trace(T)/3.0
return sph if not tensor else np.eye(3)*sph
class TestTensor:
n = 1000
c = np.random.randint(n)
@pytest.mark.parametrize('vectorized,single',[(tensor.deviatoric, deviatoric),
(tensor.eigenvalues, eigenvalues),
(tensor.eigenvectors, eigenvectors),
(tensor.symmetric, symmetric),
(tensor.transpose, transpose),
(tensor.spherical, spherical),
])
def test_vectorize_1_arg(self,vectorized,single):
epsilon = np.random.rand(self.n,3,3)
epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3))
for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)):
assert np.allclose(single(epsilon[i]),v)
def test_symmetric(self):
"""Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
x = np.random.rand(self.n,3,3)
assert np.allclose(tensor.symmetric(x)*2.0,tensor.transpose(x)+x)
def test_transpose(self):
"""Ensure that a symmetric tensor equals its transpose."""
x = tensor.symmetric(np.random.rand(self.n,3,3))
assert np.allclose(tensor.transpose(x),x)
def test_eigenvalues(self):
"""Ensure that the characteristic polynomial can be solved."""
A = tensor.symmetric(np.random.rand(self.n,3,3))
lambd = tensor.eigenvalues(A)
s = np.random.randint(self.n)
for i in range(3):
assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0)
def test_eigenvalues_and_vectors(self):
"""Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial."""
A = tensor.symmetric(np.random.rand(self.n,3,3))
lambd = tensor.eigenvalues(A)
x = tensor.eigenvectors(A)
s = np.random.randint(self.n)
for i in range(3):
assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0)
def test_eigenvectors_RHS(self):
"""Ensure that RHS coordinate system does only change sign of determinant."""
A = tensor.symmetric(np.random.rand(self.n,3,3))
LRHS = np.linalg.det(tensor.eigenvectors(A,RHS=False))
RHS = np.linalg.det(tensor.eigenvectors(A,RHS=True))
assert np.allclose(np.abs(LRHS),RHS)
def test_spherical_deviatoric_part(self):
"""Ensure that full tensor is sum of spherical and deviatoric part."""
x = np.random.rand(self.n,3,3)
assert np.allclose(tensor.spherical(x,True) + tensor.deviatoric(x),
x)
def test_spherical_mapping(self):
"""Ensure that mapping to tensor is correct."""
x = np.random.rand(self.n,3,3)
tnsr = tensor.spherical(x,True)
scalar = tensor.spherical(x,False)
assert np.allclose(np.linalg.det(tnsr),
scalar**3.0)
def test_deviatoric(self):
I_n = np.broadcast_to(np.eye(3),(self.n,3,3))
r = np.logical_not(I_n)*np.random.rand(self.n,3,3)
assert np.allclose(tensor.deviatoric(I_n+r),r)