From 9bc05acba27d8543b658fd4deec854c7c79d42cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Jan 2016 12:25:26 +0100 Subject: [PATCH] re-introduced special eigenvalues routine for 3x3 matrices --- code/math.f90 | 64 ++++++++++++++++++++++++++++++++++++-- code/plastic_dislotwin.f90 | 11 +++---- 2 files changed, 65 insertions(+), 10 deletions(-) diff --git a/code/math.f90 b/code/math.f90 index 6814d279f..a5909788b 100644 --- a/code/math.f90 +++ b/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)) & diff --git a/code/plastic_dislotwin.f90 b/code/plastic_dislotwin.f90 index 2ebd39b84..9f8b1dbc7 100644 --- a/code/plastic_dislotwin.f90 +++ b/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) = eigValues + 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, 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)