!***********************************************************************!
! TITLE: CFVR8005 - Cam-clay material
! SUBTITLE: Triaxial test of a Cam-clay material. Axisymmetric case
!
! DESCRIPTION: Conventional triaxial test of a cylindrical specimen
! DESCRIPTION:
! DESCRIPTION: References:
! DESCRIPTION:   [1] "Inelastic Analysis of Solids and Structures", M. Kojic y K.J. Bathe.
! DESCRIPTION:       Example 6.3.1: Conventonal triaxial test of Cam-Clay material
! DESCRIPTION:   [2] "Constitutive Laws for Engineering Materials", C.S. Desai y Siriwardane.
! DESCRIPTION:       Example 11-4: Critical state model for laboratory specimen
! DESCRIPTION:
! DESCRIPTION: Notes:
! DESCRIPTION:   - Conventional criterion of signs used in geotechnics has been followed:
! DESCRIPTION:            compressions and shortenings => positive
! DESCRIPTION:            tractions and stretchings    => negative
! DESCRIPTION:   - Target values for kinematic quantities are computed from the numerically
! DESCRIPTION:     obtained using different formulas to those used for the numerical ones
!
! ELEMENT TYPE: PLANE182
! MODULES:
! UNITS: SI
! KEYWORD1: Materials
! KEYWORD2: Geotechnics
! KEYWORD3: Material Models
!
!***********************************************************************!

  FINISH
  ~CFCLEAR,,1
  NomFile='CFVR8005'
  /TITLE, %NomFile%, Cam-clay material

! -------------------------------------------------------------------------------------
! Initial data
! -------------------------------------------------------------------------------------
  ~UNITS,SI
  /PREP7

  M      = 1.0      ! slope of the critical state line (CSL) in q-log(p) plane
  LAM    = 0.14     ! slope of the isotropic consolidation line
  KAP    = 0.026    ! slope of the unloading-reloading line
  VOID0  = 1.08     ! initial void ratio
  P0     = 114.E+3  ! initial preconsolidation ratio
  E0     = 9900.E+3 ! initial Young's modulus
  NU     = 0.3      ! initial Poisson's ratio
  SIG0   = 100.E+3  ! initial stresses
  SIGYMX = 2.4*SIG0 ! maximum value of the axial stress
  NTSTP  = 150      ! number of substeps

! -------------------------------------------------------------------------------------
! Model definition
! -------------------------------------------------------------------------------------
  ! Materials

  ~CFMP,1,LIB,SOIL ,,CL
  ~CFMP,1,SOIL,KPLA,,3
  ~CFMP,1,SOIL,KP0,,1
  ~CFMP,1,SOIL,M,,M
  ~CFMP,1,SOIL,LAM,,LAM
  ~CFMP,1,SOIL,KAP,,KAP
  TMP = 1. + VOID0 + LAM*LOG(P0)+ KAP*LOG(SIG0/P0)
  ~CFMP,1,SOIL,VICL,,TMP
  ~CFMP,1,SOIL,P0,,P0
  ~CFMP,1,SOIL,EXst,,E0
  ~CFMP,1,SOIL,NUXYst,,NU

  ! Element types

  ET,1,PLANE182
  KEYOPT,1,3,1

  ! Meshing

  N,1,0.,0.
  N,2,1.,0.
  N,3,1.,1.
  N,4,0.,1.
  E,1,2,3,4

  ! Loads

  D,1,ALL
  D,2,UY
  D,4,UX

  *DIM,FS1,TABLE,2,,,TIME
  *DIM,FS2,TABLE,2,,,TIME
  FS1(1,0) = 0.,1.
  FS1(1,1) = SIG0,SIG0
  FS2(1,0) = 0.,1.
  FS2(1,1) = SIG0,SIGYMX
  SFE,1,2,PRES,,%FS1%
  SFE,1,3,PRES,,%FS2%

  INISTATE,SET,DTYP,STRE
  INISTATE,DEFINE,1,,,,-SIG0,-SIG0,-SIG0,0.,0.,0.

!--------------------------------------------------------------------------------------
! DATA CHECK
!--------------------------------------------------------------------------------------
! Data comparison number
  NCOMP    = 5
  NCOMP_CH = 0

! Array dimensions
  *DIM,LABEL,CHAR,NCOMP,1
  *DIM,LABEL_CH,CHAR,NCOMP_CH,1
  *DIM,VALUE,,NCOMP,3
  *DIM,VALUE_CH,CHAR,NCOMP_CH,3
  *DIM,TOLER,,NCOMP,2

! Labels
!--------------------------------------------------------------------------------------
  LABEL(1) = 'SX'
  LABEL(2) = 'SY'
  LABEL(3) = 'EPVEL' ! elastic volumetric strain
  LABEL(4) = 'EPVPL' ! plastic volumetric strain
  LABEL(5) = 'VOL'   ! specific volume

 ! Resolution
 !--------------------------------------------------------------------------------------
  /SOLU

  NLGEOM,OFF
  AUTOTS,OFF
  NROPT,UNSYM
  LNSRCH,ON
  !NEQIT,50 Se comenta pues en la version 2019r1 se da una ostia ansys cuando se especifica el numero de iteraciones y existen variables de estado
  CUTCONTROL,PLSLIMIT,5.

  ERESX,NO
  OUTRES,NSOL,1
  OUTRES,ESOL,1
  OUTRES,SVAR,1

  TIME,1.
  NSUBST,NTSTP
  SOLVE

  /POST26

  NUMVAR,25

  ESOL,2,1,1,S,Y
  ADD,2,2,,,,,,-1.
  ESOL,3,1,1,S,X
  ADD,3,3,,,,,,-1.
  ESOL,4,1,1,EPEL,Y
  ESOL,5,1,1,EPPL,Y
  ESOL,6,1,1,EPEL,X
  ESOL,7,1,1,EPPL,X
  ADD,10,4,6,,,,,-1.,-2.
  ADD,11,5,7,,,,,-1.,-2.
  ADD,12,2,3,,,,,1/3,2/3
  ESOL,13,1,1,SVAR,1
  ESOL,14,1,1,SVAR,6

  VGET,WRK1,10
  VGET,WRK2,12
  VGET,WRK3,14
  *VOPER,WRK4,KAP,DIV,WRK2
  *VOPER,WRK4,WRK4,DIV,WRK3
  *VOPER,WRK4,WRK4,INT1,WRK2,WRK1(1)
  VPUT,WRK4,15

  VGET,WRK1,11
  VGET,WRK2,13
  VGET,WRK3,14
  TMP = LAM-KAP
  *VOPER,WRK4,TMP,DIV,WRK2
  *VOPER,WRK4,WRK4,DIV,WRK3
  *VOPER,WRK4,WRK4,INT1,WRK2,WRK1(1)
  VPUT,WRK4,16

  VGET,WRK1,13
  TMP = WRK1(1)
  *VOPER,WRK1,WRK1,DIV,TMP
  *VFUN,WRK1,LOG,WRK1
  TMP = -LAM
  *VOPER,WRK1,TMP,MULT,WRK1
  VGET,WRK2,14
  TMP = WRK2(1)
  *VOPER,WRK1,TMP,ADD,WRK1
  VPUT,WRK1,17

! Correct values
!--------------------------------------------------------------------------------------
  VALUE(1,1) = SIG0
  VALUE(2,1) = SIGYMX
  VGET,WRK,15
  VALUE(3,1) = WRK(NTSTP)
  VGET,WRK,16
  VALUE(4,1) = WRK(NTSTP)
  VGET,WRK,17
  VALUE(5,1) = WRK(NTSTP)

! Obtained values
!--------------------------------------------------------------------------------------
  VGET,WRK,3
  VALUE(1,2) = WRK(NTSTP)
  VGET,WRK,2
  VALUE(2,2) = WRK(NTSTP)
  VGET,WRK,10
  VALUE(3,2) = WRK(NTSTP)
  VGET,WRK,11
  VALUE(4,2) = WRK(NTSTP)
  VGET,WRK,14
  VALUE(5,2) = WRK(NTSTP)

! Warning and error tolerances
  EPS = 9.E-03
  EPS1 = 9.E-03
  TOLER(1,1) = EPS*ABS(VALUE(1,1)) $ TOLER(1,2) = TOLER(1,1)
  TOLER(2,1) = EPS*ABS(VALUE(2,1)) $ TOLER(2,2) = TOLER(2,1)
  TOLER(3,1) = EPS*ABS(VALUE(3,1)) $ TOLER(3,2) = TOLER(3,1)
  TOLER(4,1) = EPS*ABS(VALUE(4,1)) $ TOLER(4,2) = TOLER(4,1)
  TOLER(5,1) = EPS1*ABS(VALUE(5,1)) $ TOLER(5,2) = TOLER(5,1)

!--------------------------------------------------------------------------------------
! Results comparison
!--------------------------------------------------------------------------------------
!  TMP = $ WRK1 = $ WRK2 = $ WRK3 = $ WRK4 =
  COMPARA.MAC