!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              VRM    Feb. 2007                                                |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine vrm (params,J2d,J3d,compressibility,viscosity,pressure,             &
               plasticity_parameters,plasticity_type,is_plastic,flag_vrm_pb,   &
               straintot)

use constants
use definitions
use invariants

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|

!This routine computes the rescaled viscosity in the element in the case it is
!above the failure criterion

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

implicit none

type (parameters) params
double precision :: J2d,J3d,viscosity,pressure,compressibility,straintot
character (len=8) plasticity_type
double precision plasticity_parameters(9)
logical ::  is_plastic,flag_vrm_pb

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision :: yield,fail,theta,c,e2d,zeta
double precision :: sin_theta,cos_theta,sin_phi,cos_phi
double precision :: strain_soft_in,strain_soft_out,strain_soft_phi,phi,fact
double precision :: strain_soft_yield
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

flag_vrm_pb = .false.

!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
!write (*,*) 'init_e2d: ',params%init_e2d
!write (*,*) 'real e2d: ',sqrt(J2d)
!write (*,*) 'viscosity: ',viscosity
!write (*,*) 'is_plastic: ',is_plastic
!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
!write (*,*) ''

if (params%init_e2d) then
  e2d=params%e2d0
else
  e2d=sqrt(J2d)
endif

select case (trim(plasticity_type))

  case ('vM')
    yield = plasticity_parameters(1)
    ! strain softening
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_yield=plasticity_parameters(2)
      if (straintot.gt.strain_soft_in) then
        yield = (strain_soft_yield-yield)/(strain_soft_out-strain_soft_in)*(straintot-strain_soft_out) + strain_soft_yield
        yield = max(yield, strain_soft_yield)
      end if
    end if
    fail = 2.d0*viscosity*e2d-yield
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPI')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c     = plasticity_parameters(2)
    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*(3.d0-sin(phi)))
    k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0-sin(phi)))
    yield = max(c/compressibility,k+alpha*pressure)
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPII')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c     = plasticity_parameters(2)
    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*(3.d0+sin(phi)))
    k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*(3.d0+sin(phi)))
    yield = max(c/compressibility,k+alpha*pressure)
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPIII')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c     = plasticity_parameters(2)
    alpha = (2.d0*sin(phi))  /(sqrt(3.d0)*3.d0)
    k     = (6.d0*c*cos(phi))/(sqrt(3.d0)*3.d0)
    yield = max(c/compressibility,k+alpha*pressure)
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPIV')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c     = plasticity_parameters(2)
    alpha = tan(phi)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
    k     = (3.d0*c)/(sqrt(9.d0+12.d0*(tan(phi))**2.d0))
    yield = max(c/compressibility,k+alpha*pressure)
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPV') 
    alpha = plasticity_parameters(1)
    k     = plasticity_parameters(2)
    yield = max(k/compressibility,k+alpha*pressure)
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if
    
  case ('DPVI')
    ! Strain softening/hardening
    phi = plasticity_parameters(1)*pi/180.d0
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c     = plasticity_parameters(2)
    yield = max(c/compressibility,c+tan(phi)*pressure)
    fail  = 2.d0*viscosity*e2d-yield
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('DPVII')
    ! strain softening
    phi     = plasticity_parameters(1)
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c = plasticity_parameters(2)
    sin_phi=sin(phi)
    cos_phi=cos(phi)
    yield = max(plasticity_parameters(2)/compressibility,pressure*sin_phi+c*cos_phi)
    fail  = 2.d0*viscosity*e2d-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('MC')
    if (params%init_e2d) then
      theta = 0.d0
    else
      theta = lode_angle(J2d,J3d)
    endif
    sin_theta = sin(theta)
    cos_theta = cos(theta)
    ! strain softening
    phi = plasticity_parameters(1)
    strain_soft_in=plasticity_parameters(3)
    if (strain_soft_in.gt.0.d0) then
      strain_soft_out=plasticity_parameters(4)
      strain_soft_phi=plasticity_parameters(5)*pi/180.d0
      if (straintot.gt.strain_soft_in) then
        fact=(straintot-strain_soft_in)/(strain_soft_out-strain_soft_in)
        fact=min(fact,1.d0)
        phi=phi+(strain_soft_phi-phi)*fact
      endif
    endif
    c = plasticity_parameters(2)
    sin_phi = sin(phi)
    cos_phi = cos(phi)
    zeta      = cos_theta-sin_phi*sin_theta/sqrt3
    ! modified by Jean
    !fail      = -pressure*sin_phi+e2d*zeta-c*cos_phi
    fail      = pressure*sin_phi+2.d0*viscosity*e2d*zeta-c*cos_phi
    if (fail.gt.0.d0) then
      viscosity=0.5d0/e2d*(c*cos_phi-pressure*sin_phi)/zeta
      !if (viscosity<0.d0) then
      !  viscosity=0.d0
      !  !print *,viscosity,c*cos_phi,pressure*sin_phi
      !end if
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case ('Tresca')
    if (params%init_e2d) then
      theta = 0.d0
    else
      theta = lode_angle(J2d,J3d)
    endif
    yield     = plasticity_parameters(1)
    cos_theta = cos(theta)
    fail  = 2.d0*viscosity*e2d*cos_theta-yield 
    if (fail.gt.0.d0) then
      viscosity=0.5d0*yield/e2d/cos_theta
      is_plastic=.true.
    else
      is_plastic=.false.
    end if

  case default
    print *,plasticity_type
    call stop_run('error in vrm$')
    flag_vrm_pb = .true.

end select

!write (*,*) 'vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv'
!write (*,*) 'init_e2d: ',params%init_e2d
!write (*,*) 'e2d out: ',sqrt(J2d)
!write (*,*) 'viscosity out: ',viscosity
!write (*,*) 'is_plastic out: ',is_plastic
!write (*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
!write (*,*) ''

return
end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

subroutine compute_plastic_params (params,mat)
use definitions
use constants
implicit none
type (parameters) params
type (material) :: mat(0:params%nmat)

double precision c,k,alpha,phi
integer i,iproc,nproc,ierr

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

do i=0,params%nmat
  select case (trim(mat(i)%plasticity_type))
    case ('No')

    case ('vM')

    case ('DPI')

    case ('DPII')

    case ('DPIII')

    case ('DPIV')

    case ('DPV')

    case ('DPVI')

    case ('DPVII')

    case ('MC')

    case ('Tresca')

    case default
      call stop_run('error in compute_plastic_params$')

  end select
  if (iproc.eq.0) then
    write(8,*) 'material ',i,'plasticity_type ',mat(i)%plasticity_type
    write(8,*) 'plasticity_parameters',mat(i)%plasticity_parameters(:)
  end if
end do

return
end
