diff --git a/code/crystallite.f90 b/code/crystallite.f90
index 657cf6ead..fd70cf254 100644
--- a/code/crystallite.f90
+++ b/code/crystallite.f90
@@ -3975,14 +3975,12 @@ end function crystallite_integrateStress
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_orientations
  use math, only: &
-   math_pDecomposition, &
+   math_rotationalPart33, &
    math_RtoQ, &
    math_qConj
  use FEsolving, only: &
    FEsolving_execElem, &
    FEsolving_execIP
- use IO, only: &
-   IO_warning
  use material, only: &
    material_phase, &
    homogenization_Ngrains, &
@@ -4010,36 +4008,23 @@ subroutine crystallite_orientations
    neighboring_i, &              ! integration point index of my neighbor
    myPhase, &                    ! phase
    neighboringPhase
- real(pReal), dimension(3,3) :: &
-   U, &
-   R
  real(pReal), dimension(4) :: &
    orientation
- logical &
-   error
 
  ! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
 
- !$OMP PARALLEL DO PRIVATE(error,U,R,orientation)
+ !$OMP PARALLEL DO PRIVATE(orientation)
    do e = FEsolving_execElem(1),FEsolving_execElem(2)
      do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
        do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
 ! somehow this subroutine is not threadsafe, so need critical statement here; not clear, what exactly the problem is
          !$OMP CRITICAL (polarDecomp)
-         call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error)                       ! polar decomposition of Fe
+         orientation = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,g,i,e))))          ! rotational part from polar decomposition as quaternion
          !$OMP END CRITICAL (polarDecomp)
-         if (error) then
-           call IO_warning(650_pInt, e, i, g)
-           orientation = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal]                               ! fake orientation
-         else
-           orientation = math_RtoQ(transpose(R))
-         endif
          crystallite_rotation(1:4,g,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,g,i,e), &  ! active rotation from ori0
                                                                 orientation)                               ! to current orientation (with no symmetry)
          crystallite_orientation(1:4,g,i,e) = orientation
-       enddo
-     enddo
-   enddo
+  enddo; enddo; enddo
  !$OMP END PARALLEL DO
 
 
@@ -4063,7 +4048,7 @@ subroutine crystallite_orientations
                  crystallite_disorientation(:,n,1,i,e) = &
                    lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), &
                                             crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
-                                            lattice_structure(myPhase))                              ! calculate disorientation for given symmetry
+                                            lattice_structure(myPhase))                                    ! calculate disorientation for given symmetry
                else                                                                                        ! for neighbor with different phase
                  crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal]      ! 180 degree rotation about 100 axis
                endif
@@ -4081,8 +4066,7 @@ subroutine crystallite_orientations
          call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e)
 
        endif
-     enddo
-   enddo
+   enddo; enddo
  !$OMP END PARALLEL DO
 
 end subroutine crystallite_orientations
diff --git a/code/math.f90 b/code/math.f90
index 1b0afc773..95d398fa3 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -152,7 +152,7 @@ module math
    math_sampleGaussVar, &
    math_symmetricEulers, &
    math_spectralDecompositionSym33, &
-   math_pDecomposition, &
+   math_rotationalPart33, &
    math_invariants33, &
    math_eigenvaluesSym33, &
    math_factorial, &
@@ -1936,20 +1936,22 @@ end subroutine math_spectralDecompositionSym33
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief FE = R.U
+!> @brief rotational part from polar decomposition of tensor m
 !--------------------------------------------------------------------------------------------------
-subroutine math_pDecomposition(FE,U,R,error)
+function math_rotationalPart33(m)
+ use IO, only: &
+   IO_warning
 
  implicit none
- real(pReal), intent(in), dimension(3,3) :: FE
- real(pReal), intent(out), dimension(3,3) :: R, U
- logical, intent(out) :: error
+ real(pReal), intent(in), dimension(3,3) :: m
+ real(pReal), dimension(3,3) :: math_rotationalPart33
+ real(pReal), dimension(3,3) :: U, mSquared , Uinv, EB
  real(pReal), dimension(3) :: EV
- real(pReal), dimension(3,3) :: ce, Uinv, EB
+ logical :: error
 
- ce = math_mul33x33(math_transpose33(FE),FE)
 
- call math_spectralDecompositionSym33(ce,EV,EB,error)
+ mSquared = math_mul33x33(math_transpose33(m),m)
+ call math_spectralDecompositionSym33(mSquared,EV,EB,error)
 
  U = sqrt(EV(1)) * math_tensorproduct33(EB(1:3,1),EB(1:3,1)) &
    + sqrt(EV(2)) * math_tensorproduct33(EB(1:3,2),EB(1:3,2)) &
@@ -1957,14 +1959,13 @@ subroutine math_pDecomposition(FE,U,R,error)
 
  Uinv = math_inv33(U)
  if (all(abs(Uinv) <= tiny(Uinv))) then                                                             ! math_inv33 returns zero when failed, avoid floating point equality comparison
-   R = 0.0_pReal
-   error = .True.
+   math_rotationalPart33 = math_I3
+   call IO_warning(650_pInt)
  else
-   R = math_mul33x33(FE,Uinv)
-   error = .False. .and. error
+   math_rotationalPart33 = math_mul33x33(m,Uinv)
  endif
 
-end subroutine math_pDecomposition
+end function math_rotationalPart33
 
 
 !--------------------------------------------------------------------------------------------------