module_invariants.f90

module invariants

implicit none

contains

!this function computes the trace of a tensor
!Tr[\sigma]=\sum_{i} \sigma_{ii}

! Original version pre 05/11
!function trace (exx,eyy,ezz,exy,eyz,ezx)
!double precision trace
!double precision exx,eyy,ezz,exy,eyz,ezx
!trace=exx+eyy+ezz
!end function trace

function trace (exx,eyy,ezz)
double precision trace
double precision exx,eyy,ezz
trace=exx+eyy+ezz
end function trace

!this subroutine takes the 6 components of a tensor
!computes the trace and outputs the deviatoric tensor
!\sigma^d=\sigma-\frac{1}{3}Tr[\sigma] {\bf 1}

subroutine deviatoric (invariants_2d,exx,eyy,ezz,exy,eyz,ezx)
implicit none
logical :: invariants_2d
double precision exx,eyy,ezz,exy,eyz,ezx
double precision trace
trace=(exx+eyy+ezz)/3.d0
exx=exx-trace
eyy=eyy-trace
ezz=ezz-trace
if (invariants_2d) then
  eyy=0.d0
  exy=0.d0
  eyz=0.d0
endif
end subroutine deviatoric

!this function computes the second invariant:
!J_2=\frac{1}{2} \sum_{ij} e_{ij}e_{ji}
!in the case of a symmetric tensor

function second_invariant (exx,eyy,ezz,exy,eyz,ezx)

double precision second_invariant 
double precision exx,eyy,ezz,exy,eyz,ezx
second_invariant=(exx**2+eyy**2+ezz**2)/2.d0+exy**2+eyz**2+ezx**2
end function second_invariant


!this function computes the second invariant:
!J_3=\frac{1}{3} \sum_{ijk} e_{ij}e_{jk}e_{ji}
!in the case of a symmetric tensor

function third_invariant (exx,eyy,ezz,exy,eyz,ezx)
double precision third_invariant
double precision exx,eyy,ezz,exy,eyz,ezx
third_invariant=(exx*(     exx**2 + 3.d0*exy**2 + 3.d0*ezx**2)       &
                +eyy*(3.d0*exy**2 +      eyy**2 + 3.d0*eyz**2)       &
                +ezz*(3.d0*ezx**2 + 3.d0*eyz**2 +      ezz**2))/3.d0 &
                +2.d0*exy*eyz*ezx
end function third_invariant


!this function computes the Lode angle
!\theta=-\frac{1}{3}\sin^{-1} \left( \frac{3\sqrt{3}}{2}  J_3 / J2^{3/2}  \right)

function lode_angle (J2d,J3d)
use constants
double precision lode_angle,J2d,J3d, arg
if ((j2d.eq.0.d0).or.(j3d.eq.0.d0)) then
  lode_angle = 0.d0
else
  arg = -1.5d0*sqrt3*J3d/J2d**(1.5d0)
  if (abs(arg).gt.1.d0) then
    ! write(*,*) 'lode_angle; arg is ', arg, ', setting it to 1.d0'
    arg = 1.d0
  endif
  lode_angle=asin(arg)/3.d0
  ! lode_angle=asin(-1.5d0*sqrt3*J3d/J2d**(1.5d0))/3.d0
endif
end function lode_angle

end module invariants