diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 5595a9693..d09a1fb9d 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -956,9 +956,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 diff --git a/python/damask/_result.py b/python/damask/_result.py index 99eed6d00..06598e13b 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -944,7 +944,7 @@ class Result: @staticmethod def _add_rotational(F): return { - 'data': mechanics.rotational(F['data']), + 'data': mechanics.rotational(F['data']).as_matrix(), 'label': f"R({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index fc713401b..5fb8109c3 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -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 'As 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_Euler_angles(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__ diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 817d5aa9c..c401a3032 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,6 +1,7 @@ """Finite-strain continuum mechanics.""" from . import tensor as _tensor +from . import _rotation import numpy as _np @@ -107,11 +108,11 @@ def rotational(T): Returns ------- - R : numpy.ndarray of shape (...,3,3) - Rotational part. + R : damask.Rotation of shape (...) + Rotational part of the vector. """ - return _polar_decomposition(T,'R')[0] + return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0]) def strain(F,t,m): @@ -260,7 +261,7 @@ def _polar_decomposition(T,requested): output.append(_np.einsum('...ji,...jk',R,T)) if len(output) == 0: - raise ValueError('Output needs to be out of V,R,U') + raise ValueError('output needs to be out of V, R, U') return tuple(output) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index dea66e1ef..ce34c48de 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -265,7 +265,7 @@ class TestResult: default.add_rotational('F') loc = {'F': default.get_dataset_location('F'), 'R(F)': default.get_dataset_location('R(F)')} - in_memory = mechanics.rotational(default.read_dataset(loc['F'],0)) + in_memory = mechanics.rotational(default.read_dataset(loc['F'],0)).as_matrix() in_file = default.read_dataset(loc['R(F)'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 352d7cf37..8b99f34ec 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -682,6 +682,9 @@ 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): @@ -831,6 +834,13 @@ 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_vector, 'as_Rodrigues_vector'), diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 94ab042ed..b4672dbbf 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -3,6 +3,7 @@ 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) @@ -96,7 +97,6 @@ class TestMechanics: @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.rotational, rotational), (mechanics.stretch_left, stretch_left), (mechanics.stretch_right, stretch_right), ]) @@ -106,6 +106,14 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): assert np.allclose(single(epsilon[i]),v) + 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.rotational(epsilon_vec).as_matrix(), + mechanics.rotational(epsilon).as_matrix().shape)): + assert np.allclose(rotational(epsilon[i]),v) + + @pytest.mark.parametrize('vectorized,single',[(mechanics.stress_Cauchy, stress_Cauchy), (mechanics.stress_second_Piola_Kirchhoff, stress_second_Piola_Kirchhoff) ]) @@ -138,7 +146,7 @@ class TestMechanics: 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(F) + R = mechanics.rotational(F).as_matrix() V = mechanics.stretch_left(F) U = mechanics.stretch_right(F) assert np.allclose(np.matmul(R,U), @@ -162,7 +170,7 @@ class TestMechanics: @pytest.mark.parametrize('t',['V','U']) def test_strain_rotation(self,m,t): """Ensure that pure rotation results in no strain.""" - F = mechanics.rotational(np.random.rand(self.n,3,3)) + F = Rotation.from_random(self.n).as_matrix() assert np.allclose(mechanics.strain(F,t,m), 0.0) @@ -173,7 +181,7 @@ 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(x))), + assert np.allclose(np.abs(np.linalg.det(mechanics._polar_decomposition(x,'R')[0])), 1.0) def test_deviatoric_Mises(self): @@ -203,3 +211,7 @@ class TestMechanics: """Ensure that sherical stress has max shear of 0.0.""" 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')