Elastoplastic.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_elastoplastic
  use hecmw_util
  use mmaterial
  use m_elasticlinear
  use muyield
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
 
  implicit none
 
  private
  public :: calelastoplasticmatrix
  public :: backwardeuler
  public :: updateepstate
  public :: update_planestressj2
  public :: planestressj2tangent
  public :: checkplanestressj2material, planestressj2errormessage
  integer, parameter, public :: j2_unsupported=1, j2_bad_state=2, j2_bad_property=3
  integer, parameter, public :: j2_singular=4, j2_no_convergence=5, j2_nonfinite=6, j2_temperature=7
 
  real(kind=kreal), parameter :: id(6,6) = reshape( &
    & (/  2.d0/3.d0, -1.d0/3.d0, -1.d0/3.d0,  0.d0,  0.d0,  0.d0,   &
    &    -1.d0/3.d0,  2.d0/3.d0, -1.d0/3.d0,  0.d0,  0.d0,  0.d0,   &
    &    -1.d0/3.d0, -1.d0/3.d0,  2.d0/3.d0,  0.d0,  0.d0,  0.d0,   &
    &          0.d0,       0.d0,       0.d0, 0.5d0,  0.d0,  0.d0,   &
    &          0.d0,       0.d0,       0.d0,  0.d0, 0.5d0,  0.d0,   &
    &          0.d0,       0.d0,       0.d0,  0.d0,  0.d0, 0.5d0/), &
    & (/6, 6/))
  real(kind=kreal), parameter :: i2(6) = (/ 1.d0, 1.d0, 1.d0, 0.d0, 0.d0, 0.d0 /)
 
  integer, parameter :: vm_elastic = 0
  integer, parameter :: vm_plastic = 1
 
  integer, parameter :: mc_elastic       = 0
  integer, parameter :: mc_plastic_surf  = 1
  integer, parameter :: mc_plastic_right = 2
  integer, parameter :: mc_plastic_left  = 3
  integer, parameter :: mc_plastic_apex  = 4
 
  integer, parameter :: dp_elastic      = 0
  integer, parameter :: dp_plastic_surf = 1
  integer, parameter :: dp_plastic_apex = 2
 
  real(kind=kreal), parameter :: shell_shear_correction = 5.0d0/6.0d0
  integer, parameter :: plane_stress_components(3) = (/ 1, 2, 4 /)
 
contains
 
  logical function isplanestressj2perfectplasticity( matl )
    type(tmaterial), intent(in) :: matl
 
    isplanestressj2perfectplasticity = .false.
    if( .not. iselastoplastic(matl%mtype) ) return
    if( getelastictype(matl%mtype) /= 0 ) return
    if( getyieldfunction(matl%mtype) /= 0 ) return
    if( gethardentype(matl%mtype) /= 0 ) return
    if( abs(matl%variables(m_plconst2)) > tiny(1.0d0) ) return
    isplanestressj2perfectplasticity = .true.
  end function isplanestressj2perfectplasticity
 
  subroutine checkplanestressj2material(matl, ierr)
    type(tmaterial), intent(in) :: matl
    integer(kind=kint), intent(out) :: ierr
    type(ttable), pointer :: table
    real(kind=kreal) :: constants(4), elastic(2)
    logical :: missing
 
    ierr = j2_unsupported
    if( .not. isplanestressj2perfectplasticity(matl) ) return
    if( associated(matl%shell_var) ) then
      if( any(matl%shell_var%ortho /= 0) ) return
    endif
    call fetch_table(mc_isoelastic, matl%dict, table, missing)
    if( .not. missing ) then
      if( table%ndepends /= 0 ) then
        ierr = j2_temperature
        return
      endif
    endif
    call fetch_tabledata(mc_isoelastic, matl%dict, elastic, missing)
    ierr = j2_bad_property
    if( missing ) return
    constants = (/ elastic, matl%variables(m_plconst1), matl%variables(m_plconst2) /)
    if( .not. all(ieee_is_finite(constants)) ) return
    if( constants(1) <= 0.0d0 .or. constants(2) <= -1.0d0 .or. constants(2) >= 0.5d0 ) return
    if( constants(3) <= 0.0d0 ) return
    ierr = 0
  end subroutine checkplanestressj2material
 
  function planestressj2errormessage(ierr) result(message)
    integer(kind=kint), intent(in) :: ierr
    character(len=128) :: message
 
    select case(ierr)
    case(j2_unsupported)
      message = 'Plane-stress J2 requires isotropic elasticity and perfect plasticity'
    case(j2_bad_state)
      message = 'Invalid plane-stress J2 history state'
    case(j2_bad_property)
      message = 'Invalid or missing plane-stress J2 properties: require E > 0, -1 < nu < 0.5, yield stress > 0'
    case(j2_singular)
      message = 'Singular thickness tangent in plane-stress J2 condensation'
    case(j2_no_convergence)
      message = 'Plane-stress J2 local iteration did not converge'
    case(j2_nonfinite)
      message = 'Non-finite value in plane-stress J2 material update'
    case(j2_temperature)
      message = 'Temperature-dependent elasticity is not supported by plane-stress J2 shells'
    case default
      message = 'Unknown plane-stress J2 error'
    end select
  end function planestressj2errormessage
 
  subroutine condenseplanestressj2tangent( tangent3d, shear_stiffness, tangent, ierr )
    real(kind=kreal), intent(in) :: tangent3d(6,6)
    real(kind=kreal), intent(in) :: shear_stiffness
    real(kind=kreal), intent(out) :: tangent(5,5)
    integer(kind=kint), intent(out) :: ierr
 
    integer :: i, j, ii, jj
 
    ierr = 0
    tangent = 0.0d0
    if( .not. all(ieee_is_finite(tangent3d)) ) then
      ierr = j2_nonfinite
      return
    endif
    if( abs(tangent3d(3,3)) <= 100.0d0*epsilon(1.0d0)*maxval(abs(tangent3d)) ) then
      ierr = j2_singular
      return
    endif
 
    do i = 1, 3
      ii = plane_stress_components(i)
      do j = 1, 3
        jj = plane_stress_components(j)
        tangent(i,j) = tangent3d(ii,jj)-tangent3d(ii,3)*tangent3d(3,jj)/tangent3d(3,3)
      enddo
    enddo
    tangent(4,4) = shear_stiffness
    tangent(5,5) = shear_stiffness
  end subroutine condenseplanestressj2tangent
 
  subroutine update_planestressj2( matl, strain_increment, stress_bak, plstrain, &
      fstatus, stress, strain_increment_out, tangent, istat, fstatus_out, plpotential, energy_increment, ierr )
    type(tmaterial), intent(in) :: matl
    real(kind=kreal), intent(in) :: strain_increment(6), stress_bak(6)
    real(kind=kreal), intent(in) :: plstrain, fstatus(:)
    real(kind=kreal), intent(out) :: stress(6), strain_increment_out(6), tangent(5,5)
    integer, intent(out) :: istat
    real(kind=kreal), intent(out) :: fstatus_out(:), plpotential, energy_increment
    integer(kind=kint), intent(out) :: ierr
 
    integer, parameter :: maxiter = 25
    integer :: iter
    real(kind=kreal) :: elastic(6,6), tangent3d(6,6)
    real(kind=kreal) :: strain_increment_work(6)
    real(kind=kreal) :: stress_work(6), stress_bak_work(6)
    real(kind=kreal) :: shear_stiffness, residual, tolerance
 
    ierr = 0
    stress = 0.0d0
    strain_increment_out = strain_increment
    tangent = 0.0d0
    istat = 0
    fstatus_out = 0.0d0
    plpotential = 0.0d0
    energy_increment = 0.0d0
 
    call checkplanestressj2material(matl, ierr)
    if( ierr /= 0 ) return
    if( size(fstatus) /= size(fstatus_out) .or. size(fstatus) < 1 ) then
      ierr = j2_bad_state
      return
    endif
    if( .not. all(ieee_is_finite(strain_increment)) .or. .not. all(ieee_is_finite(stress_bak)) .or. &
        .not. all(ieee_is_finite(fstatus)) .or. .not. ieee_is_finite(plstrain) ) then
      ierr = j2_nonfinite
      return
    endif
    fstatus_out = fstatus
 
    ! Only the in-plane response enters the three-dimensional J2 update.
    strain_increment_work = strain_increment
    strain_increment_work(5:6) = 0.0d0
    stress_bak_work = stress_bak
    stress_bak_work(5:6) = 0.0d0
    strain_increment_work(3) = 0.0d0
    call calelasticmatrix( matl, d3, elastic, 0.0d0 )
    shear_stiffness = shell_shear_correction*elastic(5,5)
 
    do iter = 1, maxiter
      ! Each plane-stress iteration starts from the committed J2 state.
      fstatus_out = fstatus
      fstatus_out(1) = plstrain
      stress_work = stress_bak_work+matmul(elastic, strain_increment_work)
      call backwardeuler( matl, stress_work, plstrain, istat, fstatus_out, plpotential, 0.0d0 )
      call calelastoplasticmatrix( matl, d3, stress_work, istat, fstatus_out, plstrain, tangent3d, 0.0d0 )
 
      if( .not. all(ieee_is_finite(stress_work)) .or. .not. all(ieee_is_finite(tangent3d)) .or. &
          .not. all(ieee_is_finite(fstatus_out)) ) then
        ierr = j2_nonfinite
        return
      endif
 
      residual = stress_work(3)
      tolerance = 1.0d-10*max(matl%variables(m_plconst1), maxval(abs(stress_work)))
      if( abs(residual) <= tolerance ) exit
      if( abs(tangent3d(3,3)) <= 100.0d0*epsilon(1.0d0)*maxval(abs(tangent3d)) ) then
        ierr = j2_singular
        return
      endif
      strain_increment_work(3) = strain_increment_work(3)-residual/tangent3d(3,3)
    enddo
    if( iter > maxiter ) then
      ierr = j2_no_convergence
      return
    endif
 
    call condenseplanestressj2tangent(tangent3d, shear_stiffness, tangent, ierr)
    if( ierr /= 0 ) return
 
    stress = stress_work
    stress(3) = 0.0d0
    stress(5:6) = stress_bak(5:6)+shear_stiffness*strain_increment(5:6)
    strain_increment_out = strain_increment_work
    strain_increment_out(5:6) = strain_increment(5:6)
    ! Use the elastic predictor work before applying the plastic correction.
    energy_increment = dot_product(stress_bak_work, strain_increment_work) &
      +0.5d0*dot_product(matmul(elastic, strain_increment_work), strain_increment_work)+plpotential &
      +dot_product(stress_bak(5:6)+0.5d0*shear_stiffness*strain_increment(5:6), strain_increment(5:6))
    if( .not. ieee_is_finite(energy_increment) ) ierr = j2_nonfinite
  end subroutine update_planestressj2
 
  subroutine planestressj2tangent( matl, stress, istat, fstatus, plstrain, tangent, ierr )
    type(tmaterial), intent(in) :: matl
    real(kind=kreal), intent(in) :: stress(6), fstatus(:), plstrain
    integer, intent(in) :: istat
    real(kind=kreal), intent(out) :: tangent(5,5)
    integer(kind=kint), intent(out) :: ierr
 
    real(kind=kreal) :: tangent3d(6,6), stress_work(6), elastic(6,6)
 
    ierr = 0
    tangent = 0.0d0
    call checkplanestressj2material(matl, ierr)
    if( ierr /= 0 ) return
    if( size(fstatus) < 1 ) then
      ierr = j2_bad_state
      return
    endif
    if( .not. all(ieee_is_finite(stress)) .or. .not. all(ieee_is_finite(fstatus)) .or. &
        .not. ieee_is_finite(plstrain) ) then
      ierr = j2_nonfinite
      return
    endif
    stress_work = stress
    stress_work(5:6) = 0.0d0
    call calelasticmatrix( matl, d3, elastic, 0.0d0 )
    call calelastoplasticmatrix( matl, d3, stress_work, istat, fstatus, plstrain, tangent3d, 0.0d0 )
    call condenseplanestressj2tangent(tangent3d, shell_shear_correction*elastic(5,5), tangent, ierr)
  end subroutine planestressj2tangent
 
  subroutine calelastoplasticmatrix( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
    type( tmaterial ), intent(in) :: matl
    integer, intent(in)           :: secttype
    real(kind=kreal), intent(in)  :: stress(6) 
    real(kind=kreal), intent(in)  :: extval(:) 
    real(kind=kreal), intent(in)  :: plstrain  
    integer, intent(in)           :: istat
    real(kind=kreal), intent(out) :: d(:,:)    
    real(kind=kreal), intent(in)  :: temperature
    integer(kind=kint), intent(in), optional :: hdflag
 
    integer :: ytype,hdflag_in
 
    hdflag_in = 0
    if( present(hdflag) ) hdflag_in = hdflag
 
    ytype = getyieldfunction( matl%mtype )
    select case (ytype)
    case (0)
      call calelastoplasticmatrix_vm( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
    case (1)
      call calelastoplasticmatrix_mc( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
    case (2)
      call calelastoplasticmatrix_dp( matl, secttype, stress, istat, extval, plstrain, d, temperature, hdflag_in )
    case (3)
      call uelastoplasticmatrix( matl%variables, stress, istat, extval, plstrain, d, temperature, hdflag_in )
    end select
  end subroutine calelastoplasticmatrix
 
  subroutine calelastoplasticmatrix_vm( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
    type( tmaterial ), intent(in) :: matl
    integer, intent(in)           :: secttype
    real(kind=kreal), intent(in)  :: stress(6) 
    real(kind=kreal), intent(in)  :: extval(:) 
    real(kind=kreal), intent(in)  :: plstrain  
    integer, intent(in)           :: istat
    real(kind=kreal), intent(out) :: d(:,:)    
    real(kind=kreal), intent(in)  :: temperature
    integer(kind=kint), intent(in) :: hdflag
 
    integer :: i,j
    logical :: kinematic
    real(kind=kreal) :: dum, a(6), g, dlambda
    real(kind=kreal) :: c1,c2,c3, back(6)
    real(kind=kreal) :: j1,j2, harden, khard, devia(6)
 
    if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
 
    call calelasticmatrix( matl, secttype, d, temperature, hdflag=hdflag )
    if( istat == vm_elastic ) return
    if( hdflag == 2 ) return
 
    harden = calhardencoeff( matl, extval(1), temperature )
 
    kinematic = iskinematicharden( matl%mtype )
    khard = 0.d0
    if( kinematic ) then
      back(1:6) = extval(2:7)
      khard = calkinematicharden( matl, extval(1) )
    endif
 
    j1 = (stress(1)+stress(2)+stress(3))
    devia(1:3) = stress(1:3)-j1/3.d0
    devia(4:6) = stress(4:6)
    if( kinematic ) devia = devia-back
    j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
      dot_product( devia(4:6), devia(4:6) )
 
    a(1:6) = devia(1:6)/sqrt(2.d0*j2)
    g = d(4,4)
    dlambda = extval(1)-plstrain
    c3 = sqrt(3.d0*j2)+3.d0*g*dlambda !trial mises stress
    c1 = 6.d0*dlambda*g*g/c3
    dum = 3.d0*g+khard+harden
    c2 = 6.d0*g*g*(dlambda/c3-1.d0/dum)
 
    do i=1,6
      do j=1,6
        d(i,j) = d(i,j) - c1*id(i,j) + c2*a(i)*a(j)
      enddo
    enddo
  end subroutine calelastoplasticmatrix_vm
 
  subroutine calelastoplasticmatrix_mc( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
    use m_utilities, only : eigen3,deriv_general_iso_tensor_func_3d
    type( tmaterial ), intent(in) :: matl
    integer, intent(in)           :: secttype
    real(kind=kreal), intent(in)  :: stress(6) 
    real(kind=kreal), intent(in)  :: extval(:) 
    real(kind=kreal), intent(in)  :: plstrain  
    integer, intent(in)           :: istat
    real(kind=kreal), intent(out) :: d(:,:)    
    real(kind=kreal), intent(in)  :: temperature
    integer(kind=kint), intent(in) :: hdflag
 
    real(kind=kreal) :: g, k, harden, r2g, r2gd3, r4gd3, r2k, youngs, poisson
    real(kind=kreal) :: phi, psi, cosphi, sinphi, cotphi, sinpsi, sphsps, r2cosphi, r4cos2phi
    real(kind=kreal) :: prnstre(3), prnprj(3,3), prnstra(3)
    integer(kind=kint) :: m1, m2, m3
    real(kind=kreal) :: c1,c2, ca1, ca2, ca3, cam, cap, cd1, cd2, cd3, cdiag, coffd
    real(kind=kreal) :: ck1, ck2, ck3
    real(kind=kreal) :: dum, da, db, dc, dd, detinv
    real(kind=kreal) :: dpsdpe(3,3)
 
    if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
 
    call calelasticmatrix( matl, secttype, d, temperature, hdflag=hdflag )
    if( istat == mc_elastic ) return
    if( hdflag == 2 ) return
 
    harden = calhardencoeff( matl, extval(1), temperature )
    g = d(4,4)
    k = d(1,1)-(4.d0/3.d0)*g
    r2g = 2.d0*g
    r2k = 2.d0*k
    r2gd3 = r2g/3.d0
    r4gd3 = 2.d0*r2gd3
    youngs = 9.d0*k*g/(3.d0*k+g)
    poisson = (3.d0*k-r2g)/(6.d0*k+r2g)
 
    phi = matl%variables(m_plconst3)
    psi = matl%variables(m_plconst4)
    sinphi = sin(phi)
    cosphi = cos(phi)
    sinpsi = sin(psi)
    sphsps = sinphi*sinpsi
    r2cosphi = 2.d0*cosphi
    r4cos2phi = r2cosphi*r2cosphi
 
    call eigen3( stress, prnstre, prnprj )
    m1 = maxloc( prnstre, 1 )
    m3 = minloc( prnstre, 1 )
    if( m1 == m3 ) then
      m1 = 1; m2 = 2; m3 = 3
    else
      m2 = 6 - (m1 + m3)
    endif
 
    c1 = 4.d0*(g*(1.d0+sphsps/3.d0)+k*sphsps)
    if( istat==mc_plastic_surf ) then
      dd= c1 + r4cos2phi*harden
      cd1 = (r2g*(1.d0+sinpsi/3.d0) + r2k*sinpsi)/dd
      cd2 = (r4gd3-r2k)*sinpsi/dd
      cd3 = (r2g*(1.d0-sinpsi/3.d0) - r2k*sinpsi)/dd
      cap = 1.d0+sinphi/3.d0
      cam = 1.d0-sinphi/3.d0
      ck1 = 1.d0-2.d0*cd1*sinphi
      ck2 = 1.d0+2.d0*cd2*sinphi
      ck3 = 1.d0+2.d0*cd3*sinphi
      dpsdpe(m1,m1) = r2g*( 2.d0/3.d0-cd1*cap)+k*ck1
      dpsdpe(m1,m2) = (k-r2gd3)*ck1
      dpsdpe(m1,m3) = r2g*(-1.d0/3.d0+cd1*cam)+k*ck1
      dpsdpe(m2,m1) = r2g*(-1.d0/3.d0+cd2*cap)+k*ck2
      dpsdpe(m2,m2) = r4gd3*( 1.d0-cd2*sinphi)+k*ck2
      dpsdpe(m2,m3) = r2g*(-1.d0/3.d0-cd2*cam)+k*ck2
      dpsdpe(m3,m1) = r2g*(-1.d0/3.d0+cd3*cap)+k*ck3
      dpsdpe(m3,m2) = (k-r2gd3)*ck3
      dpsdpe(m3,m3) = r2g*( 2.d0/3.d0-cd3*cam)+k*ck3
    else if( istat==mc_plastic_apex ) then
      cotphi = cosphi/sinphi
      dpsdpe(:,:) = k*(1.d0-(k/(k+harden*cotphi*cosphi/sinpsi)))
    else ! EDGE
      if( istat==mc_plastic_right ) then
        c2 = r2g*(1.d0+sinphi+sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
      else if( istat==mc_plastic_left ) then
        c2 = r2g*(1.d0-sinphi-sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
      endif
      dum = r4cos2phi*harden
      da = c1 + dum
      db = c2 + dum
      dc = db
      dd = da
      detinv = 1.d0/(da*dd-db*dc)
      ca1 = r2g*(1.d0+sinpsi/3.d0)+r2k*sinpsi
      ca2 = (r4gd3-r2k)*sinpsi
      ca3 = r2g*(1.d0-sinpsi/3.d0)-r2k*sinpsi
      cdiag = k+r4gd3
      coffd = k-r2gd3
      if( istat==mc_plastic_right ) then
        dpsdpe(m1,m1) = cdiag+ca1*(db-dd-da+dc)*(r2g+(r2k+r2gd3)*sinphi)*detinv
        dpsdpe(m1,m2) = coffd+ca1*(r2g*(da-db)+((db-dd-da+dc)*(r2k+r2gd3)+(dd-dc)*r2g)*sinphi)*detinv
        dpsdpe(m1,m3) = coffd+ca1*(r2g*(dd-dc)+((db-dd-da+dc)*(r2k+r2gd3)+(da-db)*r2g)*sinphi)*detinv
        dpsdpe(m2,m1) = coffd+(ca2*(dd-db)+ca3*(da-dc))*(r2g+(r2k+r2gd3)*sinphi)*detinv
        dpsdpe(m2,m2) = cdiag+(ca2*((r2k*(dd-db)-(db*r2gd3+dd*r4gd3))*sinphi+db*r2g) &
            &                 +ca3*((r2k*(da-dc)+(da*r2gd3+dc*r4gd3))*sinphi-da*r2g))*detinv
        dpsdpe(m2,m3) = coffd+(ca2*((r2k*(dd-db)+(db*r4gd3+dd*r2gd3))*sinphi-dd*r2g) &
            &                 +ca3*((r2k*(da-dc)-(da*r4gd3+dc*r2gd3))*sinphi+dc*r2g))*detinv
        dpsdpe(m3,m1) = coffd+(ca2*(da-dc)+ca3*(dd-db))*(r2g+(r2k+r2gd3)*sinphi)*detinv
        dpsdpe(m3,m2) = coffd+(ca2*((r2k*(da-dc)+(da*r2gd3+dc*r4gd3))*sinphi-da*r2g) &
            &                 +ca3*((r2k*(dd-db)-(db*r2gd3+dd*r4gd3))*sinphi+db*r2g))*detinv
        dpsdpe(m3,m3) = cdiag+(ca2*((r2k*(da-dc)-(da*r4gd3+dc*r2gd3))*sinphi+dc*r2g) &
            &                 +ca3*((r2k*(dd-db)+(db*r4gd3+dd*r2gd3))*sinphi-dd*r2g))*detinv
      else if( istat==mc_plastic_left ) then
        dpsdpe(m1,m1) = cdiag+(ca1*((r2k*(db-dd)-(db*r4gd3+dd*r2gd3))*sinphi-dd*r2g) &
            &                 +ca2*((r2k*(da-dc)-(da*r4gd3+dc*r2gd3))*sinphi-dc*r2g))*detinv
        dpsdpe(m1,m2) = coffd+(ca1*((r2k*(db-dd)+(db*r2gd3+dd*r4gd3))*sinphi+db*r2g) &
            &                 +ca2*((r2k*(da-dc)+(da*r2gd3+dc*r4gd3))*sinphi+da*r2g))*detinv
        dpsdpe(m1,m3) = coffd+(ca1*(db-dd)+ca2*(da-dc))*(-r2g+(r2k+r2gd3)*sinphi)*detinv
        dpsdpe(m2,m1) = coffd+(ca1*((r2k*(dc-da)+(da*r4gd3+dc*r2gd3))*sinphi+dc*r2g) &
            &                 +ca2*((r2k*(dd-db)+(db*r4gd3+dd*r2gd3))*sinphi+dd*r2g))*detinv
        dpsdpe(m2,m2) = cdiag+(ca1*((r2k*(dc-da)-(da*r2gd3+dc*r4gd3))*sinphi-da*r2g) &
            &                 +ca2*((r2k*(dd-db)-(db*r2gd3+dd*r4gd3))*sinphi-db*r2g))*detinv
        dpsdpe(m2,m3) = coffd+(ca1*(dc-da)+ca2*(dd-db))*(-r2g+(r2k+r2gd3)*sinphi)*detinv
        dpsdpe(m3,m1) = coffd+ca3*((r2k*(-db+dd+da-dc)+(db-da)*r4gd3+(dd-dc)*r2gd3)*sinphi+(dd-dc)*r2g)*detinv
        dpsdpe(m3,m2) = coffd+ca3*((r2k*(-db+dd+da-dc)+(da-db)*r2gd3+(dc-dd)*r4gd3)*sinphi+(da-db)*r2g)*detinv
        dpsdpe(m3,m3) = cdiag+ca3*(-db+dd+da-dc)*(-r2g+(r2k+r2gd3)*sinphi)*detinv
      endif
    endif
    ! compute principal elastic strain from principal stress
    prnstra(1) = (prnstre(1)-poisson*(prnstre(2)+prnstre(3)))/youngs
    prnstra(2) = (prnstre(2)-poisson*(prnstre(1)+prnstre(3)))/youngs
    prnstra(3) = (prnstre(3)-poisson*(prnstre(1)+prnstre(2)))/youngs
    call deriv_general_iso_tensor_func_3d(dpsdpe, d, prnprj, prnstra, prnstre)
  end subroutine calelastoplasticmatrix_mc
 
  subroutine calelastoplasticmatrix_dp( matl, sectType, stress, istat, extval, plstrain, D, temperature, hdflag )
    type( tmaterial ), intent(in) :: matl
    integer, intent(in)           :: secttype
    real(kind=kreal), intent(in)  :: stress(6) 
    real(kind=kreal), intent(in)  :: extval(:) 
    real(kind=kreal), intent(in)  :: plstrain  
    integer, intent(in)           :: istat
    real(kind=kreal), intent(out) :: d(:,:)    
    real(kind=kreal), intent(in)  :: temperature
    integer(kind=kint), intent(in) :: hdflag
 
    integer :: i,j
    real(kind=kreal) :: dum, a(6), dlambda, g, k
    real(kind=kreal) :: j1,j2, eta, xi, etabar, harden, devia(6)
    real(kind=kreal) :: alpha, beta, c1, c2, c3, c4, ca, devia_norm
 
    if( secttype /=d3 ) stop "Elastoplastic calculation support only Solid element currently"
 
    call calelasticmatrix( matl, secttype, d, temperature, hdflag=hdflag )
    if( istat == dp_elastic ) return   ! elastic state
    if( hdflag == 2 ) return
 
    harden = calhardencoeff( matl, extval(1), temperature )
 
    g = d(4,4)
    k = d(1,1)-(4.d0/3.d0)*g
 
    eta = matl%variables(m_plconst3)
    xi = matl%variables(m_plconst4)
    etabar = matl%variables(m_plconst5)
 
    if( istat==dp_plastic_surf ) then
      j1 = (stress(1)+stress(2)+stress(3))
      devia(1:3) = stress(1:3)-j1/3.d0
      devia(4:6) = stress(4:6)
      j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
          dot_product( devia(4:6), devia(4:6) )
 
      devia_norm = sqrt(2.d0*j2)
      a(1:6) = devia(1:6)/devia_norm
      dlambda = extval(1)-plstrain
      ca = 1.d0 / (g + k*eta*etabar + xi*xi*harden)
      dum = sqrt(2.d0)*devia_norm
      c1 = 4.d0*g*g*dlambda/dum
      c2 = 2.d0*g*(2.d0*g*dlambda/dum - g*ca)
      c3 = sqrt(2.d0)*g*ca*k
      c4 = k*k*eta*etabar*ca
      do j=1,6
        do i=1,6
          d(i,j) = d(i,j) - c1*id(i,j) + c2*a(i)*a(j) &
              - c3*(eta*a(i)*i2(j) + etabar*i2(i)*a(j)) &
              - c4*i2(i)*i2(j)
        enddo
      enddo
    else ! istat==DP_PLASTIC_APEX
      alpha = xi/etabar
      beta = xi/eta
      c1 = k*(1.d0 - k/(k + alpha*beta*harden))
      do j=1,6
        do i=1,6
          d(i,j) = c1*i2(i)*i2(j)
        enddo
      enddo
    endif
  end subroutine calelastoplasticmatrix_dp
 
  real(kind=kreal) function calhardencoeff( matl, pstrain, temp )
    type( tmaterial ), intent(in)          :: matl
    real(kind=kreal), intent(in)           :: pstrain 
    real(kind=kreal), intent(in)           :: temp 
 
    integer :: htype
    logical :: ierr
    real(kind=kreal) :: s0, s1,s2, ef, ina(2)
 
    calhardencoeff = -1.d0
    htype = gethardentype( matl%mtype )
    select case (htype)
      case (0)  ! Linear hardening
        calhardencoeff = matl%variables(m_plconst2)
      case (1)  ! Multilinear approximation
        ina(1) = temp;  ina(2)=pstrain
        call fetch_tablegrad( mc_yield, ina, matl%dict, calhardencoeff, ierr )
      case (2)  ! Swift
        s0= matl%variables(m_plconst1)
        s1= matl%variables(m_plconst2)
        s2= matl%variables(m_plconst3)
        calhardencoeff = s1*s2*( s0+pstrain )**(s2-1)
      case (3)  ! Ramberg-Osgood
        s0= matl%variables(m_plconst1)
        s1= matl%variables(m_plconst2)
        s2= matl%variables(m_plconst3)
        ef = calcurryield( matl, pstrain, temp )
        calhardencoeff = s1*(ef/s1)**(1.d0-s2) /(s0*s2)
      case(4)   ! Prager
        calhardencoeff = 0.d0
      case(5)   ! Prager+linear
        calhardencoeff = matl%variables(m_plconst2)
    end select
  end function
 
  real(kind=kreal) function calkinematicharden( matl, pstrain )
    type( tmaterial ), intent(in) :: matl
    real(kind=kreal), intent(in)  :: pstrain 
 
    integer :: htype
    htype = gethardentype( matl%mtype )
    select case (htype)
      case(4, 5)   ! Prager
        calkinematicharden = matl%variables(m_plconst3)
      case default
        calkinematicharden = 0.d0
    end select
  end function
 
  real(kind=kreal) function calcurrkinematic( matl, pstrain )
    type( tmaterial ), intent(in) :: matl
    real(kind=kreal), intent(in)  :: pstrain 
 
    integer :: htype
    htype = gethardentype( matl%mtype )
    select case (htype)
      case(4, 5)   ! Prager
        calcurrkinematic = matl%variables(m_plconst3)*pstrain
      case default
        calcurrkinematic = 0.d0
    end select
  end function
 
  real(kind=kreal) function calcurryield( matl, pstrain, temp )
    type( tmaterial ), intent(in) :: matl
    real(kind=kreal), intent(in)  :: pstrain 
    real(kind=kreal), intent(in)  :: temp  
 
    integer :: htype
    real(kind=kreal) :: s0, s1,s2, ina(2), outa(1)
    logical :: ierr
    calcurryield = -1.d0
    htype = gethardentype( matl%mtype )
 
    select case (htype)
      case (0, 5)  ! Linear hardening, Linear+Parger hardening
        calcurryield = matl%variables(m_plconst1)+matl%variables(m_plconst2)*pstrain
      case (1)  ! Multilinear approximation
        ina(1) = temp;  ina(2)=pstrain
        call fetch_tabledata(mc_yield, matl%dict, outa, ierr, ina)
        if( ierr ) stop "Fail to get yield stress!"
        calcurryield = outa(1)
      case (2)  ! Swift
        s0= matl%variables(m_plconst1)
        s1= matl%variables(m_plconst2)
        s2= matl%variables(m_plconst3)
        calcurryield = s1*( s0+pstrain )**s2
      case (3)  ! Ramberg-Osgood
        s0= matl%variables(m_plconst1)
        s1= matl%variables(m_plconst2)
        s2= matl%variables(m_plconst3)
        if( pstrain<=s0 ) then
          calcurryield = s1
        else
          calcurryield = s1*( pstrain/s0 )**(1.d0/s2)
        endif
      case (4)  ! Parger hardening
        calcurryield = matl%variables(m_plconst1)
    end select
  end function
 
  subroutine backwardeuler( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag )
    type( tmaterial ), intent(in)    :: matl
    real(kind=kreal), intent(inout)  :: stress(6)   
    real(kind=kreal), intent(in)     :: plstrain    
    integer, intent(inout)           :: istat
    real(kind=kreal), intent(inout)  :: fstat(:)    
    real(kind=kreal), intent(inout)  :: plpotential    
    real(kind=kreal), intent(in)     :: temp  
    integer(kind=kint), intent(in), optional :: hdflag
 
    integer :: ytype, hdflag_in
 
    hdflag_in = 0
    if( present(hdflag) ) hdflag_in = hdflag
 
    ytype = getyieldfunction( matl%mtype )
    ! The return mappings leave the trial state untouched when the trial stress is elastic.
    if( ytype /= 3 ) then
      fstat(1) = plstrain
      if( iskinematicharden( matl%mtype ) ) fstat(2:7) = fstat(8:13)
    endif
    select case (ytype)
    case (0)
      call backwardeuler_vm( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag_in )
    case (1)
      call backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag_in )
    case (2)
      call backwardeuler_dp( matl, stress, plstrain, istat, fstat, temp, hdflag_in )
    case (3)
      call ubackwardeuler( matl%variables, stress, plstrain, istat, fstat, temp, hdflag_in )
    end select
  end subroutine backwardeuler
 
  subroutine backwardeuler_vm( matl, stress, plstrain, istat, fstat, plpotential, temp, hdflag )
    type( tmaterial ), intent(in)    :: matl
    real(kind=kreal), intent(inout)  :: stress(6)   
    real(kind=kreal), intent(in)     :: plstrain    
    integer, intent(inout)           :: istat
    real(kind=kreal), intent(inout)  :: fstat(:)    
    real(kind=kreal), intent(inout)  :: plpotential    
    real(kind=kreal), intent(in)     :: temp  
    integer(kind=kint), intent(in)   :: hdflag
 
    real(kind=kreal), parameter :: tol =1.d-8
    integer, parameter          :: maxiter = 10
    real(kind=kreal) :: dlambda, f
    integer :: i
    real(kind=kreal) :: youngs, poisson, pstrain, ina(1), ee(2)
    real(kind=kreal) :: j1, j2, h, kh, kk, dd, eqvs, yd, g, k, devia(6)
    logical          :: kinematic, ierr
    real(kind=kreal) :: betan, back(6)
 
    plpotential = 0.d0
    kinematic = iskinematicharden( matl%mtype )
    if( kinematic ) back(1:6) = fstat(8:13)
 
    j1 = (stress(1)+stress(2)+stress(3))
    devia(1:3) = stress(1:3)-j1/3.d0
    devia(4:6) = stress(4:6)
    if( kinematic ) devia = devia-back
    j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
      dot_product( devia(4:6), devia(4:6) )
 
    eqvs = dsqrt( 3.d0*j2 )
    yd = calcurryield( matl, plstrain, temp )
    f = eqvs - yd
 
    if( abs(f/yd)<tol ) then  ! yielded
      istat = vm_plastic
      return
    elseif( f<0.d0 ) then   ! not yielded or unloading
      istat = vm_elastic
      return
    endif
    if( hdflag == 2 ) return
 
    istat = vm_plastic      ! yielded
    kh = 0.d0; kk=0.d0; betan=0.d0
    if( kinematic ) then
      betan = calcurrkinematic( matl, plstrain )  ! keep back = alpha_n (loaded above) so it accumulates at fstat update and is restored into stress
    else
      back(:)=0.d0
    endif
 
    ina(1) = temp
    call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
    if( ierr ) then
      stop " fail to fetch young's modulus in elastoplastic calculation"
    else
      youngs = ee(1)
      poisson = ee(2)
    endif
    if( youngs==0.d0 ) stop "YOUNG's ratio==0"
    g = youngs/ ( 2.d0*(1.d0+poisson) )
    k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
 
    dlambda = 0.d0
    pstrain = plstrain
 
    do i=1,maxiter
      h= calhardencoeff( matl, pstrain, temp )
      if( kinematic ) then
        kh = calkinematicharden( matl, pstrain )
      endif
      dd= 3.d0*g+h+kh
      dlambda = dlambda+f/dd
      if( dlambda<0.d0 ) then
        dlambda = 0.d0
        pstrain = plstrain
        istat=vm_elastic; exit
      endif
      pstrain = plstrain+dlambda
      yd = calcurryield( matl, pstrain, temp )
      if( kinematic ) then
        kk = calcurrkinematic( matl, pstrain )
      endif
      f = eqvs-3.d0*g*dlambda-yd -(kk-betan)
      if( abs(f/yd)<tol ) exit
      ! if( i==MAXITER ) then
      !   stop 'ERROR: BackwardEuler_VM: convergence failure'
      ! endif
    enddo
    if( kinematic ) then
      kk = calcurrkinematic( matl, pstrain )
      fstat(2:7) = back(:)+(kk-betan)*devia(:)/eqvs
    endif
    devia(:) = (1.d0-3.d0*dlambda*g/eqvs)*devia(:)
    stress(1:3) = devia(1:3)+j1/3.d0
    stress(4:6) = devia(4:6)
    stress(:)= stress(:)+back(:)
 
    fstat(1) = pstrain
 
    h= calhardencoeff( matl, pstrain, temp ) !a
    yd = calcurryield( matl, plstrain, temp ) !b
    plpotential = -0.5d0*(eqvs-yd)*(eqvs-yd)/(h+3.d0*g)
 
  end subroutine backwardeuler_vm
 
  subroutine backwardeuler_mc( matl, stress, plstrain, istat, fstat, temp, hdflag )
    use m_utilities, only : eigen3
    type( tmaterial ), intent(in)    :: matl
    real(kind=kreal), intent(inout)  :: stress(6)   
    real(kind=kreal), intent(in)     :: plstrain    
    integer, intent(inout)           :: istat
    real(kind=kreal), intent(inout)  :: fstat(:)    
    real(kind=kreal), intent(in)     :: temp  
    integer(kind=kint), intent(in)   :: hdflag
 
    real(kind=kreal), parameter :: tol =1.d-8
    integer, parameter          :: maxiter = 10
    real(kind=kreal) :: dlambda, f, mat(3,3)
    integer :: i, m1, m2, m3
    real(kind=kreal) :: youngs, poisson, pstrain, ina(1), ee(2)
    real(kind=kreal) :: h, dd, eqvs, cohe, g, k
    real(kind=kreal) :: prnstre(3), prnprj(3,3), tstre(3,3)
    real(kind=kreal) :: phi, psi, trialprn(3)
    logical          :: ierr
    real(kind=kreal) :: c1, c2, cs1, cs2, cs3
    real(kind=kreal) :: sinphi, cosphi, sinpsi, sphsps, r2cosphi, r4cos2phi, cotphi
    real(kind=kreal) :: da, db, dc, depv, detinv, dlambdb, dum, eps, eqvsb, fb
    real(kind=kreal) :: pt, p, resid
 
    phi = matl%variables(m_plconst3)
    psi = matl%variables(m_plconst4)
    sinphi = sin(phi)
    cosphi = cos(phi)
    r2cosphi = 2.d0*cosphi
 
    call eigen3( stress, prnstre, prnprj )
    trialprn = prnstre
    m1 = maxloc( prnstre, 1 )
    m3 = minloc( prnstre, 1 )
    if( m1 == m3 ) then
      m1 = 1; m2 = 2; m3 = 3
    else
      m2 = 6 - (m1 + m3)
    endif
 
    eqvs = prnstre(m1)-prnstre(m3) + (prnstre(m1)+prnstre(m3))*sinphi
    cohe = calcurryield( matl, plstrain, temp )
    f = eqvs - r2cosphi*cohe
 
    if( abs(f/cohe)<tol ) then  ! yielded
      istat = mc_plastic_surf
      return
    elseif( f<0.d0 ) then   ! not yielded or unloading
      istat = mc_elastic
      return
    endif
    if( hdflag == 2 ) return
 
    istat = mc_plastic_surf   ! yielded
 
    ina(1) = temp
    call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
    if( ierr ) then
      stop " fail to fetch young's modulus in elastoplastic calculation"
    else
      youngs = ee(1)
      poisson = ee(2)
    endif
    if( youngs==0.d0 ) stop "YOUNG's ratio==0"
    g = youngs/ ( 2.d0*(1.d0+poisson) )
    k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
 
    dlambda = 0.d0
    pstrain = plstrain
 
    sinpsi = sin(psi)
    sphsps = sinphi*sinpsi
    r4cos2phi = r2cosphi*r2cosphi
    c1 = 4.d0*(g*(1.d0+sphsps/3.d0)+k*sphsps)
    do i=1,maxiter
      h= calhardencoeff( matl, pstrain, temp )
      dd= c1 + r4cos2phi*h
      dlambda = dlambda+f/dd
      if( r2cosphi*dlambda<0.d0 ) then
        if( cosphi==0.d0 ) stop "Math error in return mapping"
        dlambda = 0.d0
        pstrain = plstrain
        istat = mc_elastic; exit
      endif
      pstrain = plstrain + r2cosphi*dlambda
      cohe = calcurryield( matl, pstrain, temp )
      f = eqvs - c1*dlambda - r2cosphi*cohe
      if( abs(f/cohe)<tol ) exit
      ! if( i==MAXITER ) then
      !   stop 'ERROR: BackwardEuler_MC: convergence failure'
      ! endif
    enddo
    cs1 =2.d0*g*(1.d0+sinpsi/3.d0) + 2.d0*k*sinpsi
    cs2 =(4.d0*g/3.d0-2.d0*k)*sinpsi
    cs3 =2.d0*g*(1.d0-sinpsi/3.d0) - 2.d0*k*sinpsi
    prnstre(m1) = prnstre(m1)-cs1*dlambda
    prnstre(m2) = prnstre(m2)+cs2*dlambda
    prnstre(m3) = prnstre(m3)+cs3*dlambda
    eps = (abs(prnstre(m1))+abs(prnstre(m2))+abs(prnstre(m3)))*tol
    if( prnstre(m1) < prnstre(m2)-eps .or. prnstre(m2) < prnstre(m3)-eps ) then
      ! return mapping to EDGE
      prnstre = trialprn
      dlambda = 0.d0
      dlambdb = 0.d0
      if( (1.d0-sinpsi)*prnstre(m1) - 2*prnstre(m2) + (1.d0+sinpsi)*prnstre(m3) > 0) then
        istat = mc_plastic_right
        eqvsb = prnstre(m1)-prnstre(m2) + (prnstre(m1)+prnstre(m2))*sinphi
        c2 = 2.d0*g*(1.d0+sinphi+sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
      else
        istat = mc_plastic_left
        eqvsb = prnstre(m2)-prnstre(m3) + (prnstre(m2)+prnstre(m3))*sinphi
        c2 = 2.d0*g*(1.d0-sinphi-sinpsi-sphsps/3.d0) + 4.d0*k*sphsps
      endif
      cohe = calcurryield( matl, plstrain, temp )
      f = eqvs - r2cosphi*cohe
      fb = eqvsb - r2cosphi*cohe
      pstrain = plstrain
      do i=1,maxiter
        h= calhardencoeff( matl, pstrain, temp )
        dum = r4cos2phi*h
        da = c1 + dum
        db = c2 + dum
        dc = db
        dd = da
        detinv = 1.d0/(da*dd-db*dc)
        dlambda = dlambda + detinv*( dd*f - db*fb)
        dlambdb = dlambdb + detinv*(-dc*f + da*fb)
        pstrain = plstrain + r2cosphi*(dlambda+dlambdb)
        cohe = calcurryield( matl, pstrain, temp )
        f = eqvs - c1*dlambda - c2*dlambdb - r2cosphi*cohe
        fb = eqvsb - c2*dlambda - c1*dlambdb - r2cosphi*cohe
        if( (abs(f)+abs(fb))/(abs(eqvs)+abs(eqvsb)) < tol ) exit
        ! if( i==MAXITER ) then
        !   stop 'ERROR: BackwardEuler_MC: convergence failure(2)'
        ! endif
      enddo
      if( istat==mc_plastic_right ) then
        prnstre(m1) = prnstre(m1)-cs1*(dlambda+dlambdb)
        prnstre(m2) = prnstre(m2)+cs2*dlambda+cs3*dlambdb
        prnstre(m3) = prnstre(m3)+cs3*dlambda+cs2*dlambdb
      else
        prnstre(m1) = prnstre(m1)-cs1*dlambda+cs2*dlambdb
        prnstre(m2) = prnstre(m2)+cs2*dlambda-cs1*dlambdb
        prnstre(m3) = prnstre(m3)+cs3*(dlambda+dlambdb)
      endif
      eps = (abs(prnstre(m1))+abs(prnstre(m2))+abs(prnstre(m3)))*tol
      if( prnstre(m1) < prnstre(m2)-eps .or. prnstre(m2) < prnstre(m3)-eps ) then
        ! return mapping to APEX
        prnstre = trialprn
        istat = mc_plastic_apex
        if( sinphi==0.d0 ) stop 'ERROR: BackwardEuler_MC: phi==0.0'
        if( sinpsi==0.d0 ) stop 'ERROR: BackwardEuler_MC: psi==0.0'
        depv = 0.d0
        cohe = calcurryield( matl, plstrain, temp )
        cotphi = cosphi/sinphi
        pt = (stress(1)+stress(2)+stress(3))/3.d0
        resid = cotphi*cohe - pt
        pstrain = plstrain
        do i=1,maxiter
          h= calhardencoeff( matl, pstrain, temp )
          dd= cosphi*cotphi*h/sinpsi + k
          depv = depv - resid/dd
          pstrain = plstrain + cosphi*depv/sinpsi
          cohe = calcurryield( matl, pstrain,temp )
          p = pt-k*depv
          resid = cotphi*cohe-p
          if( abs(resid/cohe)<tol ) exit
          ! if( i==MAXITER ) then
          !   stop 'ERROR: BackwardEuler_MC: convergence failure(3)'
          ! endif
        enddo
        prnstre(m1) = p
        prnstre(m2) = p
        prnstre(m3) = p
      endif
    endif
    tstre(:,:) = 0.d0
    tstre(1,1)= prnstre(1); tstre(2,2)=prnstre(2); tstre(3,3)=prnstre(3)
    mat= matmul( prnprj, tstre )
    mat= matmul( mat, transpose(prnprj) )
    stress(1) = mat(1,1)
    stress(2) = mat(2,2)
    stress(3) = mat(3,3)
    stress(4) = mat(1,2)
    stress(5) = mat(2,3)
    stress(6) = mat(3,1)
 
    fstat(1) = pstrain
  end subroutine backwardeuler_mc
 
  subroutine backwardeuler_dp( matl, stress, plstrain, istat, fstat, temp, hdflag )
    type( tmaterial ), intent(in)    :: matl
    real(kind=kreal), intent(inout)  :: stress(6)   
    real(kind=kreal), intent(in)     :: plstrain    
    integer, intent(inout)           :: istat
    real(kind=kreal), intent(inout)  :: fstat(:)    
    real(kind=kreal), intent(in)     :: temp  
    integer(kind=kint), intent(in)   :: hdflag
 
    real(kind=kreal), parameter :: tol =1.d-8
    integer, parameter          :: maxiter = 10
    real(kind=kreal) :: dlambda, f
    integer :: i
    real(kind=kreal) :: youngs, poisson, pstrain, xi, ina(1), ee(2)
    real(kind=kreal) :: j1,j2,h, dd, eqvst, eqvs, cohe, g, k, devia(6), eta, etabar, pt, p
    logical          :: ierr
    real(kind=kreal) :: alpha, beta, depv, factor, resid
 
    eta = matl%variables(m_plconst3)
    xi = matl%variables(m_plconst4)
    etabar = matl%variables(m_plconst5)
 
    j1 = (stress(1)+stress(2)+stress(3))
    pt = j1/3.d0
    devia(1:3) = stress(1:3)-pt
    devia(4:6) = stress(4:6)
    j2 = 0.5d0* dot_product( devia(1:3), devia(1:3) ) +  &
      dot_product( devia(4:6), devia(4:6) )
 
    eqvst = sqrt(j2)
    cohe = calcurryield( matl, plstrain, temp )
    f = eqvst + eta*pt - xi*cohe
 
    if( abs(f/cohe)<tol ) then  ! yielded
      istat = dp_plastic_surf
      return
    elseif( f<0.d0 ) then   ! not yielded or unloading
      istat = dp_elastic
      return
    endif
    if( hdflag == 2 ) return
 
    istat = dp_plastic_surf
 
    ina(1) = temp
    call fetch_tabledata(mc_isoelastic, matl%dict, ee, ierr, ina)
    if( ierr ) then
      stop " fail to fetch young's modulus in elastoplastic calculation"
    else
      youngs = ee(1)
      poisson = ee(2)
    endif
    if( youngs==0.d0 ) stop "YOUNG's ratio==0"
    g = youngs/ ( 2.d0*(1.d0+poisson) )
    k = youngs/ ( 3.d0*(1.d0-2.d0*poisson) )
 
    dlambda = 0.d0
    pstrain = plstrain
 
    do i=1,maxiter
      h= calhardencoeff( matl, pstrain, temp )
      dd= g+k*etabar*eta+h*xi*xi
      dlambda = dlambda+f/dd
      if( xi*dlambda<0.d0 ) then
        if( xi==0.d0 ) stop "Math error in return mapping"
        dlambda = 0.d0
        pstrain = plstrain
        istat=0; exit
      endif
      pstrain = plstrain+xi*dlambda
      cohe = calcurryield( matl, pstrain, temp  )
      eqvs = eqvst-g*dlambda
      p = pt-k*etabar*dlambda
      f = eqvs + eta*p- xi*cohe
      if( abs(f/cohe)<tol ) exit
      ! if( i==MAXITER ) then
      !   stop 'ERROR: BackwardEuler_DP: convergence failure'
      ! endif
    enddo
    if( eqvs>=0.d0 ) then ! converged
      factor = 1.d0-g*dlambda/eqvst
    else                  ! return mapping to APEX
      istat = dp_plastic_apex
      if( eta==0.d0 ) stop 'ERROR: BackwardEuler_DP: eta==0.0'
      if( etabar==0.d0 ) stop 'ERROR: BackwardEuler_DP: etabar==0.0'
      alpha = xi/etabar
      beta = xi/eta
      depv=0.d0
      pstrain = plstrain
      cohe = calcurryield( matl, pstrain, temp )
      resid = beta*cohe - pt
      do i=1,maxiter
        h= calhardencoeff( matl, pstrain, temp )
        dd= alpha*beta*h + k
        depv = depv - resid/dd
        pstrain = plstrain+alpha*depv
        cohe = calcurryield( matl, pstrain, temp )
        p = pt-k*depv
        resid = beta*cohe - p
        if( abs(resid/cohe)<tol ) then
          dlambda=depv/etabar
          factor=0.d0
          exit
        endif
        ! if( i==MAXITER ) then
        !   stop 'ERROR: BackwardEuler_DP: convergence failure(2)'
        ! endif
      enddo
    endif
    devia(:) = factor*devia(:)
    stress(1:3) = devia(1:3)+p
    stress(4:6) = devia(4:6)
 
    fstat(1) = pstrain
  end subroutine backwardeuler_dp
 
  subroutine updateepstate( gauss )
    use mmechgauss
    type(tgaussstatus), intent(inout) :: gauss  ! status of curr gauss point
    gauss%plstrain= gauss%fstatus(1)
    if(iskinematicharden(gauss%pMaterial%mtype)) then
      gauss%fstatus(8:13) =gauss%fstatus(2:7)
    endif
  end subroutine
 
end module m_elastoplastic