re-introduced special eigenvalues routine for 3x3 matrices

code/math.f90

@@ -152,6 +152,7 @@ module math
    math_sampleGaussVar, &
    math_symmetricEulers, &
    math_spectralDecompositionSym33, &
+   math_spectralDecompositionSym, &
    math_rotationalPart33, &
    math_invariants33, &
    math_eigenvaluesSym33, &
@@ -1910,9 +1911,9 @@ end function math_symmetricEulers
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief eigenvalues and eigenvectors of symmetric 3x3 matrix 
+!> @brief eigenvalues and eigenvectors of symmetric matrix m
 !--------------------------------------------------------------------------------------------------
-subroutine math_spectralDecompositionSym33(M,values,vectors,error)
+subroutine math_spectralDecompositionSym(m,values,vectors,error)
 
  implicit none
  real(pReal), dimension(:,:),                  intent(in)  :: m
@@ -1931,6 +1932,63 @@ subroutine math_spectralDecompositionSym33(M,values,vectors,error)
 #endif
  error = (info == 0_pInt)
 
+end subroutine math_spectralDecompositionSym
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief eigenvalues and eigenvectors of symmetric 3x3 matrix m using an analytical expression
+!> and the general LAPACK powered version as fallback
+!> @author Joachim Kopp, Max–Planck–Institut für Kernphysik, Heidelberg (Copyright (C) 2006)
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
+!--------------------------------------------------------------------------------------------------
+subroutine math_spectralDecompositionSym33(m,values,vectors)
+ 
+ implicit none
+ 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_pInt]
+
+ T = maxval(abs(values))
+ U = MAX(T, T**2_pInt)
+ threshold = sqrt(5.0e-14_pReal * U**2_pInt)
+
+! 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_spectralDecompositionSym(m,values,vectors,error)
+   return
+ endif fallback1
+
+ 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_spectralDecompositionSym(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_crossproduct(vectors(1:3,1),vectors(1:3,2))
+
 end subroutine math_spectralDecompositionSym33
 
 
@@ -1950,7 +2008,7 @@ function math_rotationalPart33(m)
 
 
  mSquared = math_mul33x33(math_transpose33(m),m)
- call math_spectralDecompositionSym33(mSquared,EV,EB,error)
+ call math_spectralDecompositionSym33(mSquared,EV,EB)
 
  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)) &

code/plastic_dislotwin.f90

@@ -1702,7 +1702,6 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
         0, 1,-1, &
         0, 1, 1  &
         ],pReal),[ 3,6])
- logical error 
  !* Shortened notation
  of = phasememberAt(ipc,ip,el)
  ph = phaseAt(ipc,ip,el)
@@ -1783,7 +1782,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
     abs(plastic_dislotwin_sbResistance(instance)) > tiny(0.0_pReal)) then
    gdot_sb = 0.0_pReal
    dgdot_dtausb = 0.0_pReal
-   call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
+   call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors)
    do j = 1_pInt,6_pInt
      sb_s = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_sComposition(1:3,j))
      sb_m = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_mComposition(1:3,j))
@@ -2197,6 +2196,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
  use math, only: &
    pi, &
    math_Mandel6to33, &
+   math_eigenvaluesSym33, &
    math_spectralDecompositionSym33
  use material, only: &
    material_phase, &
@@ -2240,8 +2240,6 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
    gdot_slip
  real(pReal), dimension(3,3) :: eigVectors
  real(pReal), dimension (3) :: eigValues
- logical :: error
-
  
  !* Shortened notation
  of = phasememberAt(ipc,ip,el)
@@ -2517,11 +2515,10 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
          enddo ; enddo
          c = c + ns
       case (sb_eigenvalues_ID)
-        call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
+        plastic_dislotwin_postResults(c+1_pInt:c+3_pInt) = math_eigenvaluesSym33(math_Mandel6to33(Tstar_v))
-        plastic_dislotwin_postResults(c+1_pInt:c+3_pInt) = eigValues
         c = c + 3_pInt
       case (sb_eigenvectors_ID)
-        call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
+        call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors)
         plastic_dislotwin_postResults(c+1_pInt:c+9_pInt) = reshape(eigVectors,[9])
         c = c + 9_pInt
       case (stress_trans_fraction_ID)