plasticity test (phenoplus) working again with changed polar decomposition

code/math.f90

@@ -152,6 +152,8 @@ module math
    math_symmetricEulers, &
    math_spectralDecompositionSym33, &
    math_spectralDecompositionSym, &
+   math_eigenValuesVectorsSym33, &
+   math_eigenValuesVectorsSym, &
    math_rotationalPart33, &
    math_invariantsSym33, &
    math_eigenvaluesSym33, &
@@ -472,7 +474,23 @@ end function math_crossproduct
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief tensor product a \otimes b
+!> @brief tensor product A \otimes B of arbitrary sized vectors A and B
+!--------------------------------------------------------------------------------------------------
+pure function math_tensorproduct(A,B)
+
+ implicit none
+ real(pReal), dimension(:), intent(in) ::  A,B
+ real(pReal), dimension(size(A,1),size(B,1)) ::  math_tensorproduct
+
+ integer(pInt) :: i,j
+
+ forall (i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_tensorproduct(i,j) = A(i)*B(j)
+
+end function math_tensorproduct
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief tensor product A \otimes B of leght-3 vectors A and B
 !--------------------------------------------------------------------------------------------------
 pure function math_tensorproduct33(A,B)
 
@@ -682,7 +700,7 @@ pure function math_exp33(A,n)
  math_exp33 = B                                                                                     ! A^0 = eye2
 
  do i = 1_pInt,n
-   invfac = invfac/real(i)                                                                          ! invfac = 1/i!
+   invfac = invfac/real(i,pReal)                                                                    ! invfac = 1/i!
    B = math_mul33x33(B,A)
    math_exp33 = math_exp33 + invfac*B                                                               ! exp = SUM (A^i)/i!
  enddo
@@ -761,6 +779,7 @@ pure subroutine math_invert33(A, InvA, DetA, error)
  DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
 
  if (abs(DetA) <= tiny(DetA)) then
+   InvA = 0.0_pReal
    error = .true.
  else
    InvA(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
@@ -1913,7 +1932,7 @@ end function math_symmetricEulers
 !--------------------------------------------------------------------------------------------------
 !> @brief eigenvalues and eigenvectors of symmetric matrix m
 !--------------------------------------------------------------------------------------------------
-subroutine math_spectralDecompositionSym(m,values,vectors,error)
+subroutine math_eigenValuesVectorsSym(m,values,vectors,error)
 
  implicit none
  real(pReal), dimension(:,:),                  intent(in)  :: m
@@ -1924,7 +1943,7 @@ subroutine math_spectralDecompositionSym(m,values,vectors,error)
  integer(pInt) :: info
  real(pReal), dimension((64+2)*size(m,1)) :: work                                                    ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
 
- vectors = M                                                                                         ! copy matrix to input (doubles as output) array
+ vectors = m                                                                                         ! copy matrix to input (doubles as output) array
 #if(FLOAT==8)
  call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info)
 #elif(FLOAT==4)
@@ -1932,68 +1951,85 @@ subroutine math_spectralDecompositionSym(m,values,vectors,error)
 #endif
  error = (info == 0_pInt)
 
-end subroutine math_spectralDecompositionSym
+end subroutine math_eigenValuesVectorsSym
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression
+!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m
-!> and the general LAPACK powered version for arbritrary sized matrices 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)
+subroutine math_eigenValuesVectorsSym33(m,values,vectors)
- 
+
  implicit none
- real(pReal), dimension(3,3),intent(in)  :: m
+ real(pReal), dimension(3,3), intent(in)  :: m
- real(pReal), dimension(3),  intent(out) :: values
+ real(pReal), dimension(3),   intent(out) :: values
- real(pReal), dimension(3,3),intent(out) :: vectors
+ real(pReal), dimension(3,3), intent(out) :: vectors
- real(pReal) :: T, U, norm, threshold
  logical :: error
+ integer(pInt) :: info
+ real(pReal), dimension((64+2)*3) :: work                                                            ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
 
- values = math_eigenvaluesSym33(m)
+ vectors = m                                                                                         ! copy matrix to input (doubles as output) array
+#if(FLOAT==8)
+ call dsyev('V','U',3,vectors,3,values,work,(64+2)*3,info)
+#elif(FLOAT==4)
+ call ssyev('V','U',3,vectors,3,values,work,(64+2)*3,info)
+#endif
+ error = (info == 0_pInt)
 
- vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
+end subroutine math_eigenValuesVectorsSym33
-                    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
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief rotational part from polar decomposition of tensor m
+!> @brief eigenvalues and eigenvectors of symmetric matrix m
+!--------------------------------------------------------------------------------------------------
+function math_spectralDecompositionSym(m)
+
+ implicit none
+ real(pReal), dimension(:,:),      intent(in)  :: m
+ real(pReal), dimension(size(m,1))             :: values
+ real(pReal), dimension(size(m,1),size(m,1))   :: vectors
+ real(pReal), dimension(size(m,1),size(m,1))   :: math_spectralDecompositionSym
+ logical :: error
+ integer(pInt) :: i
+
+ math_spectralDecompositionSym = 0.0_pReal
+ call math_eigenValuesVectorsSym(m,values,vectors,error)
+ if(error) return
+ 
+ do i=1_pInt, size(m,1)
+   math_spectralDecompositionSym = math_spectralDecompositionSym &
+                                 +  sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i))
+ enddo
+
+end function math_spectralDecompositionSym
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m
+!--------------------------------------------------------------------------------------------------
+function math_spectralDecompositionSym33(m)
+
+ implicit none
+ real(pReal), dimension(3,3), intent(in)  :: m
+ real(pReal), dimension(3,3)              :: math_spectralDecompositionSym33
+ real(pReal), dimension(3)                :: values
+ real(pReal), dimension(3,3)              :: vectors
+ logical :: error
+ integer(pInt) :: i
+
+ math_spectralDecompositionSym33 = 0.0_pReal
+ call math_eigenValuesVectorsSym33(m,values,vectors)
+ if(error) return
+ 
+ do i=1_pInt, 3_pInt
+   math_spectralDecompositionSym33 = math_spectralDecompositionSym33 &
+                                   +  sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i))
+ enddo
+
+end function math_spectralDecompositionSym33
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief rotational part from polar decomposition of 33 tensor m
 !--------------------------------------------------------------------------------------------------
 function math_rotationalPart33(m)
  use IO, only: &
@@ -2002,17 +2038,11 @@ function math_rotationalPart33(m)
  implicit none
  real(pReal), intent(in), dimension(3,3) :: m
  real(pReal), dimension(3,3) :: math_rotationalPart33
- real(pReal), dimension(3,3) :: U, mTm , Uinv, EB
+ real(pReal), dimension(3,3) :: U , Uinv
- real(pReal), dimension(3) :: EV
-
- mTm = math_mul33x33(math_transpose33(m),m)
- call math_spectralDecompositionSym33(mTm,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)) &
-   + sqrt(EV(3)) * math_tensorproduct33(EB(1:3,3),EB(1:3,3))
 
+ U = math_spectralDecompositionSym33(math_mul33x33(transpose(m),m))
  Uinv = math_inv33(U)
+
  if (all(abs(Uinv) <= tiny(Uinv))) then                                                             ! math_inv33 returns zero when failed, avoid floating point equality comparison
    math_rotationalPart33 = math_I3
    call IO_warning(650_pInt)
@@ -2675,4 +2705,4 @@ real(pReal) pure function math_limit(a, left, right)
 
 end function math_limit
  
-end module math
+end module math

code/plastic_dislotwin.f90

@@ -1637,7 +1637,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
    math_Plain3333to99, &
    math_Mandel6to33, &
    math_Mandel33to6, &
-   math_spectralDecompositionSym, &
+   math_eigenValuesVectorsSym, &
    math_tensorproduct33, &
    math_symmetric33, &
    math_mul33x3
@@ -1783,7 +1783,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_spectralDecompositionSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
+   call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
    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,8 +2197,8 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
  use math, only: &
    pi, &
    math_Mandel6to33, &
-   math_eigenvaluesSym33, &
+   math_eigenValuesSym33, &
-   math_spectralDecompositionSym33
+   math_eigenValuesVectorsSym33
  use material, only: &
    material_phase, &
    phase_plasticityInstance,& 
@@ -2519,7 +2519,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
         plastic_dislotwin_postResults(c+1_pInt:c+3_pInt) = math_eigenvaluesSym33(math_Mandel6to33(Tstar_v))
         c = c + 3_pInt
       case (sb_eigenvectors_ID)
-        call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors)
+        call math_eigenValuesVectorsSym33(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)
@@ -2539,4 +2539,4 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
  enddo
 end function plastic_dislotwin_postResults
 
-end module plastic_dislotwin
+end module plastic_dislotwin