diff --git a/PRIVATE b/PRIVATE
index 25b7b7895..8bd09b551 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 25b7b78951ab9fb682aaf048d3e0f0d09a3be695
+Subproject commit 8bd09b5511d1e0e0ea288b47d16ce4924d75adcd
diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py
index 6a02cca08..cb5361599 100755
--- a/processing/post/addOrientations.py
+++ b/processing/post/addOrientations.py
@@ -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)
diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index 9d0350337..a1ccfa3f6 100644
--- a/python/damask/__init__.py
+++ b/python/damask/__init__.py
@@ -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
diff --git a/python/damask/_geom.py b/python/damask/_geom.py
index b4806a008..92f560e42 100644
--- a/python/damask/_geom.py
+++ b/python/damask/_geom.py
@@ -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'')
diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index f91cc9850..05301561f 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -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)
diff --git a/python/damask/_result.py b/python/damask/_result.py
index c36411ef8..4e44728da 100644
--- a/python/damask/_result.py
+++ b/python/damask/_result.py
@@ -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'],
diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index 60dd055c6..5fb8109c3 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -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
diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py
index 942be4f1f..2f4d63791 100644
--- a/python/damask/_vtk.py
+++ b/python/damask/_vtk.py
@@ -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()
diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py
index c81399d94..66f5c9916 100644
--- a/python/damask/mechanics.py
+++ b/python/damask/mechanics.py
@@ -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)))
diff --git a/python/damask/seeds.py b/python/damask/seeds.py
index 717fc620b..1c4150c2d 100644
--- a/python/damask/seeds.py
+++ b/python/damask/seeds.py
@@ -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
 
diff --git a/python/damask/tensor.py b/python/damask/tensor.py
new file mode 100644
index 000000000..00fc7e72a
--- /dev/null
+++ b/python/damask/tensor.py
@@ -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)
diff --git a/python/tests/conftest.py b/python/tests/conftest.py
index 9410ba15a..719fa8e61 100644
--- a/python/tests/conftest.py
+++ b/python/tests/conftest.py
@@ -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],
diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py
index 54ffc4742..3f9855ad4 100644
--- a/python/tests/test_Geom.py
+++ b/python/tests/test_Geom.py
@@ -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)
diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py
index ccd757e06..edc5a79d8 100644
--- a/python/tests/test_Orientation.py
+++ b/python/tests/test_Orientation.py
@@ -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,)})\
diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py
index c48f5acb8..787721339 100644
--- a/python/tests/test_Result.py
+++ b/python/tests/test_Result.py
@@ -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
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index fe90e2aea..cc22aba4d 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -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
diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py
index 48dce4707..af4ecf918 100644
--- a/python/tests/test_VTK.py
+++ b/python/tests/test_VTK.py
@@ -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)
diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py
index 5fe79de28..667747f8e 100644
--- a/python/tests/test_mechanics.py
+++ b/python/tests/test_mechanics.py
@@ -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')
diff --git a/python/tests/test_tensor.py b/python/tests/test_tensor.py
new file mode 100644
index 000000000..58194130e
--- /dev/null
+++ b/python/tests/test_tensor.py
@@ -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)