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

no need to pass determinant if not needed

This commit is contained in:
Martin Diehl 2022-09-10 10:31:18 +02:00
parent 472ef1d127
commit b11ca71786
2 changed files with 13 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/math.f90
View File

@ -486,22 +486,25 @@ end function math_inv33
!> @details Direct Cramer inversion of matrix A. Also returns determinant !> @details Direct Cramer inversion of matrix A. Also returns determinant
! Returns an error if not possible, i.e. if determinant is close to zero ! Returns an error if not possible, i.e. if determinant is close to zero
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(InvA, DetA, error, A) pure subroutine math_invert33(InvA,DetA,error, A)
real(pReal), dimension(3,3), intent(out) :: InvA real(pReal), dimension(3,3), intent(out) :: InvA
real(pReal), intent(out) :: DetA real(pReal), intent(out), optional :: DetA
logical, intent(out) :: error logical, intent(out) :: error
real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3,3), intent(in) :: A
real(pReal) :: Det
InvA(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2) InvA(1,1) = A(2,2) * A(3,3) - A(2,3) * A(3,2)
InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1) InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1)
InvA(3,1) = A(2,1) * A(3,2) - A(2,2) * A(3,1) InvA(3,1) = A(2,1) * A(3,2) - A(2,2) * A(3,1)
DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1) Det = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
if (dEq0(DetA)) then if (dEq0(Det)) then
InvA = 0.0_pReal InvA = 0.0_pReal
if (present(DetA)) DetA = 0.0_pReal
error = .true. error = .true.
else else
InvA(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2) InvA(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
@ -512,7 +515,8 @@ pure subroutine math_invert33(InvA, DetA, error, A)
InvA(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1) InvA(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1)
InvA(3,3) = A(1,1) * A(2,2) - A(1,2) * A(2,1) InvA(3,3) = A(1,1) * A(2,2) - A(1,2) * A(2,1)
InvA = InvA/DetA InvA = InvA/Det
if (present(DetA)) DetA = Det
error = .false. error = .false.
end if end if

9
src/phase_mechanical.f90
View File

@ -399,8 +399,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
real(pReal) steplengthLp, & real(pReal) steplengthLp, &
steplengthLi, & steplengthLi, &
atol_Lp, & atol_Lp, &
atol_Li, & atol_Li
devNull
integer NiterationStressLp, & ! number of stress integrations integer NiterationStressLp, & ! number of stress integrations
NiterationStressLi, & ! number of inner stress integrations NiterationStressLi, & ! number of inner stress integrations
ierr, & ! error indicator for LAPACK ierr, & ! error indicator for LAPACK
@ -417,9 +416,9 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess
Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess
call math_invert33(invFp_current,devNull,error,subFp0) call math_invert33(invFp_current,error=error,A=subFp0)
if (error) return ! error if (error) return ! error
call math_invert33(invFi_current,devNull,error,subFi0) call math_invert33(invFi_current,error=error,A=subFi0)
if (error) return ! error if (error) return ! error
A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp A = matmul(F,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
@ -541,7 +540,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
end do LiLoop end do LiLoop
invFp_new = matmul(invFp_current,B) invFp_new = matmul(invFp_current,B)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,error=error,A=invFp_new)
if (error) return ! error if (error) return ! error
phase_mechanical_P(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) phase_mechanical_P(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))