diff --git a/src/math.f90 b/src/math.f90
index 48c736f76..aa81a0ed3 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -841,14 +841,14 @@ pure function math_66toSym3333(m66,weighted)
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
     w = INVNRMMANDEL
-  endif
+  end if
 
   do i=1,6; do j=1,6
     math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
     math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j)
     math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j)
     math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j)
-  enddo; enddo
+  end do; end do
 
 end function math_66toSym3333
 
@@ -866,7 +866,6 @@ pure function math_Voigt6to33_stress(sigma_tilde) result(sigma)
                    sigma_tilde(6), sigma_tilde(2), sigma_tilde(4), &
                    sigma_tilde(5), sigma_tilde(4), sigma_tilde(3)],[3,3])
 
-
 end function math_Voigt6to33_stress
 
 
@@ -883,10 +882,40 @@ pure function math_Voigt6to33_strain(epsilon_tilde) result(epsilon)
                      0.5_pReal*epsilon_tilde(6),           epsilon_tilde(2), 0.5_pReal*epsilon_tilde(4), &
                      0.5_pReal*epsilon_tilde(5), 0.5_pReal*epsilon_tilde(4),           epsilon_tilde(3)],[3,3])
 
-
 end function math_Voigt6to33_strain
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Convert 3x3 tensor into 6 Voigt stress vector.
+!--------------------------------------------------------------------------------------------------
+pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
+
+  real(pReal), dimension(6) :: sigma_tilde
+  real(pReal), dimension(3,3), intent(in) :: sigma
+
+
+  sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), &
+                 sigma(3,2), sigma(3,1), sigma(1,2)] 
+
+end function math_33toVoigt6_stress
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Convert 3x3 tensor into 6 Voigt strain vector.
+!--------------------------------------------------------------------------------------------------
+pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
+
+  real(pReal), dimension(6) :: epsilon_tilde
+  real(pReal), dimension(3,3), intent(in) :: epsilon
+
+
+  epsilon_tilde = [          epsilon(1,1),           epsilon(2,2),           epsilon(3,3), &
+                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] 
+
+end function math_33toVoigt6_strain
+
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
 !--------------------------------------------------------------------------------------------------
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index 5792ea3cd..bb0422bd1 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -148,15 +148,17 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
     dS_dFi                                                                                          !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
 
   real(pReal), dimension(3,3) :: E
+  real(pReal), dimension(6,6) :: C66
   real(pReal), dimension(3,3,3,3) :: C
   integer :: &
     i, j
 
 
-  C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en))
+  C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)
+  C = math_Voigt66to3333(C66)
 
   E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3)                                                  !< Green-Lagrange strain in unloaded configuration
-  S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi))                                        !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
+  S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
 
   do i =1,3; do j=1,3
     dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi)))                !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko