From 5d0900ee2e7bd50316e7de4bd2ce2025f19b466a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 26 Feb 2016 15:36:24 +0100 Subject: [PATCH] plasticity test (phenoplus) working again with changed polar decomposition --- code/math.f90 | 160 ++++++++++++++++++++++--------------- code/plastic_dislotwin.f90 | 12 +-- 2 files changed, 101 insertions(+), 71 deletions(-) diff --git a/code/math.f90 b/code/math.f90 index 47c7da175..918035c9e 100644 --- a/code/math.f90 +++ b/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 -!> 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) +!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m !-------------------------------------------------------------------------------------------------- -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), intent(out) :: values - real(pReal), dimension(3,3),intent(out) :: vectors - real(pReal) :: T, U, norm, threshold + real(pReal), dimension(3,3), intent(in) :: m + real(pReal), dimension(3), intent(out) :: values + real(pReal), dimension(3,3), intent(out) :: vectors 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), & - 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 +end subroutine math_eigenValuesVectorsSym33 !-------------------------------------------------------------------------------------------------- -!> @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) :: 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)) + real(pReal), dimension(3,3) :: U , Uinv + 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 \ No newline at end of file diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index 532312bfd..e84321dee 100644 --- a/code/plastic_dislotwin.f90 +++ b/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_spectralDecompositionSym33 + math_eigenValuesSym33, & + 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 \ No newline at end of file +end module plastic_dislotwin