diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py
index 9ae1d2ed3..ee2462287 100644
--- a/python/damask/_orientation.py
+++ b/python/damask/_orientation.py
@@ -1103,7 +1103,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
@@ -1133,7 +1133,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):
@@ -1157,7 +1157,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):
@@ -1183,9 +1183,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):
@@ -1227,8 +1227,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/_rotation.py b/python/damask/_rotation.py
index e60b14ecb..4cf8cfcca 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -553,7 +553,7 @@ class Rotation:
             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)) \
diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py
index a83f27370..ef2434e00 100644
--- a/python/damask/mechanics.py
+++ b/python/damask/mechanics.py
@@ -60,7 +60,7 @@ def Cauchy(P,F):
         Cauchy stress.
 
     """
-    sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F)
+    sigma = _np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)
     return tensor.symmetric(sigma)
 
 
@@ -79,7 +79,7 @@ def deviatoric_part(T):
         Deviatoric part of T.
 
     """
-    return T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T))
+    return T - _np.einsum('...ij,...',_np.eye(3),spherical_part(T))
 
 
 def maximum_shear(T_sym):
@@ -157,7 +157,7 @@ def PK2(P,F):
         Second Piola-Kirchhoff stress.
 
     """
-    S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
+    S = _np.einsum('...jk,...kl',_np.linalg.inv(F),P)
     return tensor.symmetric(S)
 
 
@@ -199,7 +199,7 @@ def spherical_part(T,tensor=False):
 
     """
     sph = _np.trace(T,axis2=-2,axis1=-1)/3.0
-    return _np.einsum('...jk,...->...jk',_np.eye(3),sph) if tensor else sph
+    return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph
 
 
 def strain(F,t,m):
@@ -231,14 +231,12 @@ def strain(F,t,m):
         w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F))
 
     if   m > 0.0:
-        eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n))
-                                  - _np.einsum('...jk->...jk',_np.eye(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('...j,...kj->...jk',w**m,n))
-                                  + _np.einsum('...jk->...jk',_np.eye(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('...j,...kj->...jk',0.5*_np.log(w),n))
+        eps = _np.einsum('...j,...kj,...lj',0.5*_np.log(w),n,n)
 
     return eps
 
@@ -293,22 +291,22 @@ def _polar_decomposition(T,requested):
 
     """
     u, _, vh = _np.linalg.svd(T)
-    R = _np.einsum('...ij,...jk->...ik',u,vh)
+    R = _np.einsum('...ij,...jk',u,vh)
 
     output = []
     if 'R' in requested:
         output.append(R)
     if 'V' in requested:
-        output.append(_np.einsum('...ij,...kj->...ik',T,R))
+        output.append(_np.einsum('...ij,...kj',T,R))
     if 'U' in requested:
-        output.append(_np.einsum('...ji,...jk->...ik',R,T))
+        output.append(_np.einsum('...ji,...jk',R,T))
 
     return tuple(output)
 
 
 def _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
     ----------
@@ -319,7 +317,7 @@ def _Mises(T_sym,s):
 
     """
     d = deviatoric_part(T_sym)
-    return _np.sqrt(s*_np.einsum('...jk->...',d**2.0))
+    return _np.sqrt(s*_np.sum(d**2.0,axis=(-1,-2)))
 
 
 # for compatibility
diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py
index 78b441e18..0b4393ed3 100644
--- a/python/tests/test_mechanics.py
+++ b/python/tests/test_mechanics.py
@@ -127,9 +127,8 @@ class TestMechanics:
         for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)):
             assert np.allclose(single(epsilon[i]),v)
 
-    @pytest.mark.parametrize('vectorized,single',[
-        (mechanics.Cauchy,Cauchy),
-        (mechanics.PK2   ,PK2   )
+    @pytest.mark.parametrize('vectorized,single',[(mechanics.Cauchy,Cauchy),
+                                                  (mechanics.PK2   ,PK2   )
                                                  ])
     def test_vectorize_2_arg(self,vectorized,single):
         P     = np.random.rand(self.n,3,3)
@@ -227,7 +226,7 @@ class TestMechanics:
                            scalar**3.0)
 
     def test_spherical_Mises(self):
-        """Ensure that Mises equivalent strrain of spherical strain is 0."""
+        """Ensure that Mises equivalent strain 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),