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

better readable, tested, and following standard notation

This commit is contained in:
Martin Diehl 2020-03-01 23:00:06 +01:00
parent 6dfc48f89e
commit 6701af6425
2 changed files with 23 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/homogenization_mech_RGC.f90
View File

@ -669,7 +669,7 @@ module procedure mech_RGC_updateState
nDef = 0.0_pReal nDef = 0.0_pReal
do i = 1,3; do j = 1,3 do i = 1,3; do j = 1,3
do k = 1,3; do l = 1,3 do k = 1,3; do l = 1,3
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
enddo; enddo enddo; enddo
nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal ! compute the norm of the mismatch tensor
enddo; enddo enddo; enddo
@ -689,7 +689,7 @@ module procedure mech_RGC_updateState
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha & rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
*surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) & *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
*cosh(prm%ciAlpha*nDefNorm) & *cosh(prm%ciAlpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) & *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
*tanh(nDefNorm/xSmoo_RGC) *tanh(nDefNorm/xSmoo_RGC)
enddo; enddo;enddo; enddo enddo; enddo;enddo; enddo
enddo interfaceLoop enddo interfaceLoop

32
src/math.f90
View File

@ -271,31 +271,30 @@ end function math_identity4th
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief permutation tensor e_ijk used for computing cross product of two tensors !> @brief permutation tensor e_ijk
! e_ijk = 1 if even permutation of ijk ! e_ijk = 1 if even permutation of ijk
! e_ijk = -1 if odd permutation of ijk ! e_ijk = -1 if odd permutation of ijk
! e_ijk = 0 otherwise ! e_ijk = 0 otherwise
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_civita(i,j,k) real(pReal) pure function math_LeviCivita(i,j,k)
integer, intent(in) :: i,j,k integer, intent(in) :: i,j,k
math_civita = 0.0_pReal if (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
if (((i == 1).and.(j == 2).and.(k == 3)) .or. & math_LeviCivita = +1.0_pReal
((i == 2).and.(j == 3).and.(k == 1)) .or. & elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal math_LeviCivita = -1.0_pReal
if (((i == 1).and.(j == 3).and.(k == 2)) .or. & else
((i == 2).and.(j == 1).and.(k == 3)) .or. & math_LeviCivita = 0.0_pReal
((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal endif
end function math_civita end function math_LeviCivita
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief kronecker delta function d_ij !> @brief kronecker delta function d_ij
! d_ij = 1 if i = j ! d_ij = 1 if i = j
! d_ij = 0 otherwise ! d_ij = 0 otherwise
! inspired by http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function math_delta(i,j) real(pReal) pure function math_delta(i,j)
@ -1309,6 +1308,7 @@ subroutine unitTest
sort_out_ = reshape([-1,-1, +1,+5, +5,+6, +3,-2],[2,4]) sort_out_ = reshape([-1,-1, +1,+5, +5,+6, +3,-2],[2,4])
integer, dimension(5) :: range_out_ = [1,2,3,4,5] integer, dimension(5) :: range_out_ = [1,2,3,4,5]
integer, dimension(3) :: ijk
real(pReal) :: det real(pReal) :: det
real(pReal), dimension(3) :: v3_1,v3_2,v3_3,v3_4 real(pReal), dimension(3) :: v3_1,v3_2,v3_3,v3_4
@ -1425,6 +1425,16 @@ subroutine unitTest
if(math_binomial(49,6) /= 13983816) & if(math_binomial(49,6) /= 13983816) &
call IO_error(0,ext_msg='math_binomial') call IO_error(0,ext_msg='math_binomial')
ijk = cshift([1,2,3],int(r*1.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
call IO_error(0,ext_msg='math_LeviCivita(even)')
ijk = cshift([3,2,1],int(r*2.0e2_pReal))
if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
call IO_error(0,ext_msg='math_LeviCivita(odd)')
ijk = cshift([2,2,1],int(r*2.0e2_pReal))
if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3))))&
call IO_error(0,ext_msg='math_LeviCivita')
end subroutine unitTest end subroutine unitTest
end module math end module math