c
c Code modified from J-C model code supplied by ThirdWave Systems 
c     + Johnson-Cook material model
c     + Isotropic-hardening, but thermal-softening and strain rate sensitivity removed
c     + To support user-defined-yield-surface
c
c     + Material parameters:
c     UMATPAR01    A           d(1),   initial yield stress
c     UMATPAR02    B           d(2),   hardening stress coeff.
c     UMATPAR03    n           d(3),   hardening power coeff.
c     UMATPAR04    0           d(4),   no T effect
c     UMATPAR05    0           d(5),   no T effect
c     UMATPAR06    0           d(6),   no T effect
c     UMATPAR07    0           d(7),   no e rate effect
c     UMATPAR08    pre-strain  d(8),   pre-strain
c     UMATPAR09    cut-off strain d(9),cut-off to strain hardening
c
c     $Id: mat_user_yield_surface.f,v 1.8 2012/10/12 18:10:32 xhuang Exp $


      SUBROUTINE udys_hello(model_name)
!DEC$ ATTRIBUTES DLLEXPORT, ALIAS:"UDYSHELLO"::udys_hello
      implicit none
      character*(*) model_name
      model_name = "Johnson-Cook Model Version #003"
      END SUBROUTINE udys_hello

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      SUBROUTINE mat_user_yield_surface(d, yield, eps, eps_rate, temp, 
     &   dt, user_s, nvar, kStep)
c
c Sample code for user defined yield surface
c      
c  d(1:50)  UMATPAR01-UMATPAR50     Material parameters
c  yield                            Yield stress [Pa]
c  temp                             Temperature [c]
c  eps                              Equivalent plastic strain [1]
c  eps_rate                         Equibalent plastic strain rate [1/s]
c  dt                               Incremental time step [s]
c  user_s(1:nvar)                   User state variables
c  kStep                            Time step count
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

#ifndef linux
!DEC$ ATTRIBUTES DLLEXPORT :: mat_user_yield_surface
#endif

      implicit none
      double precision d(*), yield, eps, eps_rate, temp, dt, user_s(*)
      integer nvar, kStep
      double precision strain_hardening
      
      ! Initialization of user state variables
      if( kStep.eq.1 ) then
        user_s(1:nvar) = 0.d0
      endif
      
      yield = strain_hardening(d, eps)
     
!!      user_s(1) = eps       ! Save the equivalent plastic strain
     
      return 
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      SUBROUTINE mat_user_yield_surface_deriv(d, dyield, eps, eps_rate, 
     &       temp, dt, user_s, nvar, kStep)
c
c Sample code for user defined yield surface
c      
c  d(1:50)  UMATPAR01-UMATPAR50     Material parameters
c  dyield                           Yield stress [Pa]
c  temp                             Temperature [c]
c  eps                              Equivalent plastic strain [1]
c  eps_rate                         Equibalent plastic strain rate [1/s]
c  dt                               Incremental time step [s]
c  user_s(1:nvar)                   User state variables
c  kStep                            Time step count
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

#ifndef linux
!DEC$ ATTRIBUTES DLLEXPORT :: mat_user_yield_surface_deriv
#endif

      implicit none
      double precision d(*), dyield, eps, eps_rate, temp, dt, user_s(*)
      integer nvar, kStep
      double precision strain_hardening,strain_hardening_deriv
      
      dyield = strain_hardening_deriv(d, eps)
      return 
      end

      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      double precision function strain_hardening(d, alpha)
c
c Sample code for Johnson-Cook strain hardening
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      implicit none
      double precision d(*), alpha,c, alph1
      double precision A_, B_, power_n
      
      A_ = d(1)     ! initial yield stress
      B_ = d(2)     ! hardening 
      power_n = d(3)   ! power coefficient
	 c = d(9)
       alph1=alpha+d(8)
	 if(alph1.gt.c) then
 	 strain_hardening = A_ + B_ * c ** power_n
	 else
      strain_hardening = A_ + B_ * alph1 ** power_n
	 end if
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      double precision function strain_hardening_deriv(d, alpha)
c
c Sample code for calculating the derivative of the Johnson-Cook 
c strain hardening function w.r.t equivalent plastic strain
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      implicit none
      double precision d(*), alpha, c, alph1
      double precision A_, B_, power_n     
      A_ = d(1)     ! initial yield stress
      B_ = d(2)     ! hardening 
      power_n = d(3)   ! power coefficient
	 c = d(9)
       alph1=alpha+d(8)
	 if (alph1.gt.c) then
	 strain_hardening_deriv = 0.0d0
	 else
      strain_hardening_deriv = power_n * B_ * alph1 ** (power_n - 1.0d0)
	 end if
      return
      end


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      

      SUBROUTINE mat_user_yield_surface_update_state(d, eps, eps_rate,
     &  temp, dt, user_s, nvar, sigma, kStep)
c
c Sample code for user defined yield surface
c      
c  d(1:50)  UMATPAR01-UMATPAR50     Material parameters
c  yield                            Yield stress [Pa]
c  temp                             Temperature [c]
c  eps                              Equivalent plastic strain [1]
c  eps_rate                         Equibalent plastic strain rate [1/s]
c  dt                               Incremental time step [s]
c  user_s(1:nvar)                   User state variables
c  kStep                            Time step count
c  sigma(1:6)                       Stress tensor
c               sigma(1)            Sxx
c               sigma(2)            Syy
c               sigma(3)            Szz
c               sigma(4)            Sxy
c               sigma(5)            Syz
c               sigma(6)            Szx
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

#ifndef linux
!DEC$ ATTRIBUTES DLLEXPORT :: mat_user_yield_surface_update_state
#endif

      implicit none
      double precision d(*), yield, eps, eps_rate, temp, dt, user_s(*), 
     &      sigma(*) 
      integer nvar, kStep
      
      double precision pressure
      pressure = -(sigma(1) + sigma(2) + sigma(3)) / 3.d0
      user_s(2) = pressure * 1.d-6  ! pressure in MPa
      user_s(5) = sigma(1) * 1.d-6  ! Sxx in MPa
      
           
      return 
      end