diff --git a/src/math.f90 b/src/math.f90
index cfe82914f..80c3ed54e 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -72,7 +72,12 @@ module math
       3,2, &
       3,3  &
       ],[2,9])                                                                                      !< arrangement in Plain notation
-     
+
+
+  interface math_mul33xx33
+    module procedure math_tensordot
+  end interface math_mul33xx33
+
 !---------------------------------------------------------------------------------------------------
  private :: &
    unitTest
@@ -88,7 +93,7 @@ subroutine math_init
   real(pReal), dimension(4) :: randTest
   integer :: randSize
   integer, dimension(:), allocatable :: randInit
-  
+
   write(6,'(/,a)') ' <<<+-  math init  -+>>>'; flush(6)
 
   call random_seed(size=randSize)
@@ -140,7 +145,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
   else
     e = ubound(a,2)
   endif
-  
+
   if(present(sortDim)) then
     d = sortDim
   else
@@ -160,12 +165,12 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
   !> @brief Partitioning required for quicksort
   !-------------------------------------------------------------------------------------------------
   integer function qsort_partition(a, istart, iend, sort)
-  
+
     integer, dimension(:,:), intent(inout) :: a
     integer,                 intent(in)    :: istart,iend,sort
     integer, dimension(size(a,1))          :: tmp
     integer :: i,j
-   
+
     do
       ! find the first element on the right side less than or equal to the pivot point
       do j = iend, istart, -1
@@ -187,7 +192,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
         a(:,j) = tmp
       endif cross
     enddo
-  
+
   end function qsort_partition
 
 end subroutine math_sort
@@ -196,7 +201,7 @@ end subroutine math_sort
 !--------------------------------------------------------------------------------------------------
 !> @brief vector expansion
 !> @details takes a set of numbers (a,b,c,...) and corresponding multiples (x,y,z,...)
-!> to return a vector of x times a, y times b, z times c, ... 
+!> to return a vector of x times a, y times b, z times c, ...
 !--------------------------------------------------------------------------------------------------
 pure function math_expand(what,how)
 
@@ -204,9 +209,9 @@ pure function math_expand(what,how)
   integer,       dimension(:), intent(in) :: how
   real(pReal), dimension(sum(how)) ::  math_expand
   integer :: i
- 
+
   if (sum(how) == 0) return
- 
+
   do i = 1, size(how)
     math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
   enddo
@@ -346,15 +351,15 @@ end function math_inner
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 33xx33 = 1 (double contraction --> ij * ij)
+!> @brief 3x3 tensor double contraction: ij * ij
 !--------------------------------------------------------------------------------------------------
-real(pReal) pure function math_mul33xx33(A,B)
+real(pReal) pure function math_tensordot(A,B)
 
   real(pReal), dimension(3,3), intent(in) :: A,B
 
-  math_mul33xx33 = sum(A*B)
+  math_tensordot = sum(A*B)
 
-end function math_mul33xx33
+end function math_tensordot
 
 
 !--------------------------------------------------------------------------------------------------
@@ -365,7 +370,7 @@ pure function math_mul3333xx33(A,B)
   real(pReal), dimension(3,3,3,3), intent(in) :: A
   real(pReal), dimension(3,3),     intent(in) :: B
   real(pReal), dimension(3,3) :: math_mul3333xx33
-  integer :: i,j 
+  integer :: i,j
 
   do i=1,3; do j=1,3
     math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
@@ -399,7 +404,7 @@ pure function math_exp33(A,n)
   real(pReal), dimension(3,3), intent(in) :: A
   integer,                     intent(in), optional :: n
   real(pReal), dimension(3,3) :: B, math_exp33
-  
+
   real(pReal) :: invFac
   integer     :: n_,i
 
@@ -424,14 +429,14 @@ end function math_exp33
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Cramer inversion of 33 matrix (function)
-!> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e. 
+!> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
 ! if determinant is close to zero
 !--------------------------------------------------------------------------------------------------
 pure function math_inv33(A)
 
   real(pReal), dimension(3,3), intent(in) :: A
   real(pReal), dimension(3,3) :: math_inv33
-  
+
   real(pReal) :: DetA
   logical     :: error
 
@@ -526,12 +531,12 @@ subroutine math_invert(InvA, error, A)
     dgetrf, &
     dgetri
 
-  invA = A 
+  invA = A
   call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
   error = (ierr /= 0)
   call dgetri(size(A,1),InvA,size(A,1),ipiv,work,size(work,1),ierr)
   error = error .or. (ierr /= 0)
- 
+
 end subroutine math_invert
 
 
@@ -618,7 +623,7 @@ end function math_trace33
 real(pReal) pure function math_det33(m)
 
   real(pReal), dimension(3,3), intent(in) :: m
-  
+
   math_det33 = m(1,1)* (m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
              - m(1,2)* (m(2,1)*m(3,3)-m(2,3)*m(3,1)) &
              + m(1,3)* (m(2,1)*m(3,2)-m(2,2)*m(3,1))
@@ -632,7 +637,7 @@ end function math_det33
 real(pReal) pure function math_detSym33(m)
 
   real(pReal), dimension(3,3), intent(in) :: m
- 
+
   math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) &
                   + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3)
 
@@ -646,9 +651,9 @@ pure function math_33to9(m33)
 
   real(pReal), dimension(9)               :: math_33to9
   real(pReal), dimension(3,3), intent(in) :: m33
-  
+
   integer :: i
-  
+
   do i = 1, 9
     math_33to9(i) = m33(MAPPLAIN(1,i),MAPPLAIN(2,i))
   enddo
@@ -663,7 +668,7 @@ pure function math_9to33(v9)
 
   real(pReal), dimension(3,3)           :: math_9to33
   real(pReal), dimension(9), intent(in) :: v9
-  
+
   integer :: i
 
   do i = 1, 9
@@ -684,10 +689,10 @@ pure function math_sym33to6(m33,weighted)
   real(pReal), dimension(6)               :: math_sym33to6
   real(pReal), dimension(3,3), intent(in) :: m33                                                     !< symmetric matrix (no internal check)
   logical,     optional,       intent(in) :: weighted                                                !< weight according to Mandel (.true. by default)
-  
+
   real(pReal), dimension(6) :: w
   integer :: i
-  
+
   if(present(weighted)) then
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
@@ -696,7 +701,7 @@ pure function math_sym33to6(m33,weighted)
 
   do i = 1, 6
     math_sym33to6(i) = w(i)*m33(MAPNYE(1,i),MAPNYE(2,i))
-  enddo 
+  enddo
 
 end function math_sym33to6
 
@@ -715,7 +720,7 @@ pure function math_6toSym33(v6,weighted)
 
   real(pReal), dimension(6) :: w
   integer :: i
-  
+
   if(present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
@@ -778,7 +783,7 @@ pure function math_sym3333to66(m3333,weighted)
 
   real(pReal), dimension(6) :: w
   integer :: i,j
-  
+
   if(present(weighted)) then
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
@@ -803,10 +808,10 @@ pure function math_66toSym3333(m66,weighted)
   real(pReal), dimension(3,3,3,3)            :: math_66toSym3333
   real(pReal), dimension(6,6),    intent(in) :: m66
   logical,     optional,          intent(in) :: weighted                                             !< weight according to Mandel (.true. by default)
-  
+
   real(pReal), dimension(6) :: w
   integer :: i,j
-  
+
   if(present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
@@ -906,48 +911,48 @@ end subroutine math_eigenValuesVectorsSym
 ! ToDo: has wrong oder of arguments
 !--------------------------------------------------------------------------------------------------
 subroutine math_eigenValuesVectorsSym33(m,values,vectors)
- 
+
   real(pReal), dimension(3,3),intent(in)  :: m
   real(pReal), dimension(3),  intent(out) :: values
   real(pReal), dimension(3,3),intent(out) :: vectors
   real(pReal) :: T, U, norm, threshold
   logical :: error
- 
+
   values = math_eigenvaluesSym33(m)
- 
+
   vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
                      m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), &
                      m(1, 2)**2]
- 
+
   T = maxval(abs(values))
   U = max(T, T**2)
   threshold = sqrt(5.68e-14_pReal * U**2)
- 
+
  ! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2
   vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), &
                      vectors(2,2) + m(2, 3) * values(1), &
                      (m(1,1) - values(1)) * (m(2,2) - values(1)) - vectors(3,2)]
   norm = norm2(vectors(1:3, 1))
- 
+
   fallback1: if(norm < threshold) then
     call math_eigenValuesVectorsSym(m,values,vectors,error)
     return
   endif fallback1
- 
-  vectors(1:3,1) = vectors(1:3, 1) / norm 
-  
+
+  vectors(1:3,1) = vectors(1:3, 1) / norm
+
  ! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2
   vectors(1:3,2) = [ vectors(1,2) + m(1, 3) * values(2), &
                      vectors(2,2) + m(2, 3) * values(2), &
                      (m(1,1) - values(2)) * (m(2,2) - values(2)) - vectors(3,2)]
   norm = norm2(vectors(1:3, 2))
- 
+
   fallback2: if(norm < threshold) then
     call math_eigenValuesVectorsSym(m,values,vectors,error)
     return
   endif fallback2
   vectors(1:3,2) = vectors(1:3, 2) / norm
- 
+
  ! Calculate third eigenvector according to  v[2] = v[0] x v[1]
   vectors(1:3,3) = math_cross(vectors(1:3,1),vectors(1:3,2))
 
@@ -965,11 +970,11 @@ function math_eigenvectorBasisSym(m)
   real(pReal), dimension(size(m,1),size(m,1))   :: math_eigenvectorBasisSym
   logical :: error
   integer :: i
- 
+
   math_eigenvectorBasisSym = 0.0_pReal
   call math_eigenValuesVectorsSym(m,values,vectors,error)
   if(error) return
-  
+
   do i=1, size(m,1)
     math_eigenvectorBasisSym = math_eigenvectorBasisSym &
                              + sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i))
@@ -989,13 +994,13 @@ pure function math_eigenvectorBasisSym33(m)
   real(pReal) :: P, Q, rho, phi
   real(pReal), parameter :: TOL=1.e-14_pReal
   real(pReal), dimension(3,3,3) :: N, EB
- 
+
   invariants = math_invariantsSym33(m)
   EB = 0.0_pReal
- 
+
   P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
   Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
- 
+
   threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then
     values = invariants(1)/3.0_pReal
   ! this is not really correct, but at least the basis is correct
@@ -1013,7 +1018,7 @@ pure function math_eigenvectorBasisSym33(m)
     N(1:3,1:3,1) = m-values(1)*math_I3
     N(1:3,1:3,2) = m-values(2)*math_I3
     N(1:3,1:3,3) = m-values(3)*math_I3
-    twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then 
+    twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
       EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
                                                 ((values(3)-values(1))*(values(3)-values(2)))
       EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
@@ -1021,7 +1026,7 @@ pure function math_eigenvectorBasisSym33(m)
       EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
                                                 ((values(1)-values(2))*(values(1)-values(3)))
       EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
-    elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues 
+    elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
       EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
                                                 ((values(2)-values(1))*(values(2)-values(3)))
       EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
@@ -1077,7 +1082,7 @@ pure function math_eigenvectorBasisSym33_log(m)
     N(1:3,1:3,1) = m-values(1)*math_I3
     N(1:3,1:3,2) = m-values(2)*math_I3
     N(1:3,1:3,3) = m-values(3)*math_I3
-    twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then 
+    twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
       EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
                                                 ((values(3)-values(1))*(values(3)-values(2)))
       EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
@@ -1085,7 +1090,7 @@ pure function math_eigenvectorBasisSym33_log(m)
       EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
                                                 ((values(1)-values(2))*(values(1)-values(3)))
       EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
-    elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues 
+    elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
       EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
                                                 ((values(2)-values(1))*(values(2)-values(3)))
       EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
@@ -1133,7 +1138,7 @@ end function math_rotationalPart33
 ! will return NaN on error
 !--------------------------------------------------------------------------------------------------
 function math_eigenvaluesSym(m)
- 
+
   real(pReal), dimension(:,:),                  intent(in)  :: m
   real(pReal), dimension(size(m,1))                         :: math_eigenvaluesSym
   real(pReal), dimension(size(m,1),size(m,1))               :: m_
@@ -1218,7 +1223,7 @@ integer pure function math_binomial(n,k)
 
   integer, intent(in) :: n, k
   integer :: i, k_, n_
- 
+
   k_ = min(k,n-k)
   n_ = n
   math_binomial = merge(1,0,k_>-1)                                                                  ! handling special cases k < 0 or k > n
@@ -1241,7 +1246,7 @@ integer pure function math_multinomial(alpha)
   math_multinomial = 1
   do i = 1, size(alpha)
     math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
-  enddo  
+  enddo
 
 end function math_multinomial
 
@@ -1346,7 +1351,7 @@ subroutine unitTest
   call random_number(v9)
   if(any(dNeq(math_33to9(math_9to33(v9)),v9))) &
     call IO_error(0,ext_msg='math_33to9/math_9to33')
-    
+
   call random_number(t99)
   if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) &
     call IO_error(0,ext_msg='math_3333to99/math_99to3333')
@@ -1375,7 +1380,7 @@ subroutine unitTest
   call random_number(t33)
   if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
     call IO_error(0,ext_msg='math_det33/math_detSym33')
-  
+
   if(any(dNeq0(math_identity2nd(3),math_inv33(math_I3)))) &
     call IO_error(0,ext_msg='math_inv33(math_I3)')
 
@@ -1407,18 +1412,18 @@ subroutine unitTest
   if(any(dNeq0(txx_2,txx) .or. e)) &
     call IO_error(0,ext_msg='math_invert(txx)/math_identity2nd')
 
-  call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix 
+  call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix
   if(any(dNeq0(matmul(t99_2,t99)-math_identity2nd(9),tol=1.0e-9_pReal)) .or. e) &
     call IO_error(0,ext_msg='math_invert(t99)')
 
   if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
-     call IO_error(0,ext_msg='math_clip')
+    call IO_error(0,ext_msg='math_clip')
 
   if(math_factorial(10) /= 3628800) &
-     call IO_error(0,ext_msg='math_factorial')
+    call IO_error(0,ext_msg='math_factorial')
 
   if(math_binomial(49,6) /= 13983816) &
-     call IO_error(0,ext_msg='math_binomial')
+    call IO_error(0,ext_msg='math_binomial')
 
 end subroutine unitTest