simplified and got MPI Heidelberg solution for eigenvalues/vectors back
This commit is contained in:
Martin Diehl
parent
ca3e1f0da0
commit
c7ab5a9396
169
code/math.f90
169
code/math.f90
|
|
@ -150,8 +150,8 @@ module math
|
|||
math_sampleFiberOri, &
|
||||
math_sampleGaussVar, &
|
||||
math_symmetricEulers, &
|
||||
math_spectralDecompositionSym33, &
|
||||
math_spectralDecompositionSym, &
|
||||
math_eigenvectorBasisSym33, &
|
||||
math_eigenvectorBasisSym, &
|
||||
math_eigenValuesVectorsSym33, &
|
||||
math_eigenValuesVectorsSym, &
|
||||
math_rotationalPart33, &
|
||||
|
|
@ -1955,84 +1955,112 @@ end subroutine math_eigenValuesVectorsSym
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m
|
||||
!> @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)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
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), 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
|
||||
integer(pInt) :: info
|
||||
real(pReal), dimension((64+2)*3) :: 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
|
||||
#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)
|
||||
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.68e-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_eigenValuesVectorsSym(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_eigenValuesVectorsSym(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_eigenValuesVectorsSym33
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief eigenvalues and eigenvectors of symmetric matrix m
|
||||
!> @brief eigenvector basis of symmetric matrix m
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function math_spectralDecompositionSym(m)
|
||||
function math_eigenvectorBasisSym(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
|
||||
real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym
|
||||
logical :: error
|
||||
integer(pInt) :: i
|
||||
|
||||
math_spectralDecompositionSym = 0.0_pReal
|
||||
math_eigenvectorBasisSym = 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))
|
||||
math_eigenvectorBasisSym = math_eigenvectorBasisSym &
|
||||
+ sqrt(values(i)) * math_tensorproduct(vectors(:,i),vectors(:,i))
|
||||
enddo
|
||||
|
||||
end function math_spectralDecompositionSym
|
||||
end function math_eigenvectorBasisSym
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m
|
||||
!> @brief eigenvector basis of symmetric 33 matrix m
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function math_spectralDecompositionSym33(m)
|
||||
function math_eigenvectorBasisSym33(m)
|
||||
|
||||
implicit none
|
||||
real(pReal), dimension(3,3) :: math_spectralDecompositionSym33
|
||||
real(pReal), dimension(3) :: invariants, values,C
|
||||
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33
|
||||
real(pReal), dimension(3) :: invariants, values
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
real(pReal) :: P, Q, rho, phi
|
||||
real(pReal), parameter :: TOL=1.e-14_pReal
|
||||
real(pReal), dimension(3,3) :: M1, M2, M3,EB1, EB2, EB3
|
||||
real(pReal) :: D1, D2, D3
|
||||
real(pReal), dimension(3,3,3) :: N, EB
|
||||
|
||||
invariants = math_invariantsSym33(m)
|
||||
EB = 0.0_pReal
|
||||
|
||||
P=invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
|
||||
Q=-2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+invariants(1)*invariants(2)/3.0_pReal-invariants(3)
|
||||
P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
|
||||
Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
|
||||
|
||||
EB1=0.0_pReal
|
||||
EB2=0.0_pReal
|
||||
EB3=0.0_pReal
|
||||
if((ABS(P) < TOL).AND.(ABS(Q) < TOL)) then ! EV_2 = EV_1 = EV_3
|
||||
threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then
|
||||
values = invariants(1)/3.0_pReal
|
||||
! this is not really correct, but this way U is calculated
|
||||
! correctly in PDECOMPOSITION (correct is EB?=I)
|
||||
EB1(1,1)=1.0_pReal
|
||||
EB2(2,2)=1.0_pReal
|
||||
EB3(3,3)=1.0_pReal
|
||||
else
|
||||
! this is not really correct, but at least the basis is correct
|
||||
EB(1,1,1)=1.0_pReal
|
||||
EB(2,2,2)=1.0_pReal
|
||||
EB(3,3,3)=1.0_pReal
|
||||
else threeSimilarEigenvalues
|
||||
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
|
||||
phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
|
||||
values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
|
||||
|
|
@ -2040,35 +2068,36 @@ function math_spectralDecompositionSym33(m)
|
|||
cos((phi+2.0_pReal*PI)/3.0_pReal), &
|
||||
cos((phi+4.0_pReal*PI)/3.0_pReal) &
|
||||
] + invariants(1)/3.0_pReal
|
||||
C = abs([values(1)-values(2),values(2)-values(3),values(3)-values(1)])
|
||||
M1=M-values(1)*math_I3
|
||||
M2=M-values(2)*math_I3
|
||||
M3=M-values(3)*math_I3
|
||||
if (C(1) < TOL) then ! EV_2 = EV_1, no contribution from EV_2
|
||||
D3=1.0_pReal/(values(3)-values(1))/(values(3)-values(2))
|
||||
EB3=math_mul33x33(M1,M2)*D3
|
||||
EB1=math_I3-EB3
|
||||
elseif (C(2) < TOL) then ! EV_2 = EV_3, no contribution from EV_3
|
||||
D1=1.0_pReal/(values(1)-values(2))/(values(1)-values(3))
|
||||
EB1=math_mul33x33(M2,M3)*D1
|
||||
EB2=math_I3-EB1
|
||||
elseif(C(3) < TOL) then ! EV_1 = EV_3, no contribution from EV_3
|
||||
D2=1.0_pReal/(values(2)-values(1))/(values(2)-values(3))
|
||||
EB2=math_mul33x33(M1,M3)*D2
|
||||
EB1=math_I3-EB2
|
||||
else ! all three eigenvectors are different
|
||||
D1=1.0_pReal/(values(1)-values(2))/(values(1)-values(3))
|
||||
D2=1.0_pReal/(values(2)-values(1))/(values(2)-values(3))
|
||||
D3=1.0_pReal/(values(3)-values(1))/(values(3)-values(2))
|
||||
EB1=math_mul33x33(M2,M3)*D1
|
||||
EB2=math_mul33x33(M1,M3)*D2
|
||||
EB3=math_mul33x33(M1,M2)*D3
|
||||
endif
|
||||
endif
|
||||
N(1:3,1:3,1) = m-values(1)*math_I3
|
||||
N(1:3,1:3,2) = m-values(2)*math_I3
|
||||
N(1:3,1:3,3) = m-values(3)*math_I3
|
||||
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
|
||||
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
|
||||
((values(3)-values(1))*(values(3)-values(2)))
|
||||
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
|
||||
elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
|
||||
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
|
||||
((values(1)-values(2))*(values(1)-values(3)))
|
||||
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
|
||||
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
|
||||
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
|
||||
((values(2)-values(1))*(values(2)-values(3)))
|
||||
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
|
||||
else twoSimilarEigenvalues
|
||||
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
|
||||
((values(1)-values(2))*(values(1)-values(3)))
|
||||
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
|
||||
((values(2)-values(1))*(values(2)-values(3)))
|
||||
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
|
||||
((values(3)-values(1))*(values(3)-values(2)))
|
||||
endif twoSimilarEigenvalues
|
||||
endif threeSimilarEigenvalues
|
||||
|
||||
math_spectralDecompositionSym33 = sqrt(values(1)) * EB1 + sqrt(values(2)) * EB2 + sqrt(values(3)) * EB3
|
||||
math_eigenvectorBasisSym33 = sqrt(values(1)) * EB(1:3,1:3,1) &
|
||||
+ sqrt(values(2)) * EB(1:3,1:3,2) &
|
||||
+ sqrt(values(3)) * EB(1:3,1:3,3)
|
||||
|
||||
end function math_spectralDecompositionSym33
|
||||
end function math_eigenvectorBasisSym33
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -2083,7 +2112,7 @@ function math_rotationalPart33(m)
|
|||
real(pReal), dimension(3,3) :: math_rotationalPart33
|
||||
real(pReal), dimension(3,3) :: U , Uinv
|
||||
|
||||
U = math_spectralDecompositionSym33(math_mul33x33(transpose(m),m))
|
||||
U = math_eigenvectorBasisSym33(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
|
||||
|
|
@ -2143,7 +2172,7 @@ function math_eigenvaluesSym33(m)
|
|||
P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
|
||||
Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
|
||||
|
||||
if(any(abs([p,q]) < TOL)) then
|
||||
if(all(abs([P,Q]) < TOL)) then
|
||||
math_eigenvaluesSym33 = math_eigenvaluesSym(m)
|
||||
else
|
||||
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
|
||||
|
|
|
|||
|
|
@ -1207,13 +1207,11 @@ subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
|
|||
math_Plain3333to99, &
|
||||
math_Mandel6to33, &
|
||||
math_Mandel33to6, &
|
||||
math_spectralDecompositionSym33, &
|
||||
math_symmetric33, &
|
||||
math_mul33x3
|
||||
use material, only: &
|
||||
material_phase, &
|
||||
phase_plasticityInstance, &
|
||||
!plasticState, &
|
||||
phaseAt, phasememberAt
|
||||
use lattice, only: &
|
||||
lattice_Sslip, &
|
||||
|
|
|
|||
Loading…
Reference in New Issue