223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

intention more clear

This commit is contained in:
Martin Diehl 2020-03-15 14:52:37 +01:00
parent 66302fa6da
commit 65decfc48a
1 changed files with 17 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/math.f90
View File

@ -975,7 +975,6 @@ pure function math_eigenvectorBasisSym33(m)
real(pReal), dimension(3,3,3) :: N, EB
I = math_invariantsSym33(m)
EB = 0.0_pReal
P = I(2)-I(1)**2.0_pReal/3.0_pReal
Q = -2.0_pReal/27.0_pReal*I(1)**3.0_pReal+product(I(1:2))/3.0_pReal-I(3)
@ -983,6 +982,7 @@ pure function math_eigenvectorBasisSym33(m)
threeSimilarEigVals: if(all(abs([P,Q]) < TOL)) then
v = I(1)/3.0_pReal
! this is not really correct, but at least the basis is correct
EB = 0.0_pReal
EB(1,1,1)=1.0_pReal
EB(2,2,2)=1.0_pReal
EB(3,3,3)=1.0_pReal
@ -999,12 +999,15 @@ pure function math_eigenvectorBasisSym33(m)
twoSimilarEigVals: if(abs(v(1)-v(2)) < TOL) then
EB(1:3,1:3,3) = matmul(N(1:3,1:3,1),N(1:3,1:3,2))/((v(3)-v(1))*(v(3)-v(2)))
EB(1:3,1:3,1) = math_I3-EB(1:3,1:3,3)
EB(1:3,1:3,2) = 0.0_pReal
elseif (abs(v(2)-v(3)) < TOL) then twoSimilarEigVals
EB(1:3,1:3,1) = matmul(N(1:3,1:3,2),N(1:3,1:3,3))/((v(1)-v(2))*(v(1)-v(3)))
EB(1:3,1:3,2) = math_I3-EB(1:3,1:3,1)
EB(1:3,1:3,3) = 0.0_pReal
elseif (abs(v(3)-v(1)) < TOL) then twoSimilarEigVals
EB(1:3,1:3,2) = matmul(N(1:3,1:3,3),N(1:3,1:3,1))/((v(2)-v(3))*(v(2)-v(1)))
EB(1:3,1:3,1) = math_I3-EB(1:3,1:3,2)
EB(1:3,1:3,3) = 0.0_pReal
else twoSimilarEigVals
EB(1:3,1:3,1) = matmul(N(1:3,1:3,2),N(1:3,1:3,3))/((v(1)-v(2))*(v(1)-v(3)))
EB(1:3,1:3,2) = matmul(N(1:3,1:3,1),N(1:3,1:3,3))/((v(2)-v(3))*(v(2)-v(1)))
@ -1025,17 +1028,17 @@ end function math_eigenvectorBasisSym33
function math_rotationalPart(m)
real(pReal), intent(in), dimension(3,3) :: m
real(pReal), dimension(3,3) :: math_rotationalPart33
real(pReal), dimension(3,3) :: math_rotationalPart
real(pReal), dimension(3,3) :: U , Uinv
U = math_eigenvectorBasisSym33(matmul(transpose(m),m))
Uinv = math_inv33(U)
inversionFailed: if (all(dEq0(Uinv))) then
math_rotationalPart33 = math_I3
math_rotationalPart = math_I3
call IO_warning(650)
else inversionFailed
math_rotationalPart33 = matmul(m,Uinv)
math_rotationalPart = matmul(m,Uinv)
endif inversionFailed
end function math_rotationalPart