!***********************************************************************
! TITLE: CFVR8003 - Mohr-Coulomb
! SUBTITLE: CivilFEM's Mohr-Coulomb equivalent to ANSYS's Hyperbolic Extended Drucker-Prager
!
! DESCRIPTION: Checking equivalence between the material models Mohr-Coulomb of
! DESCRIPTION: CivilFEM and Hyperbolic Linear Extended Drucker-Prager of ANSYS
! DESCRIPTION:
! DESCRIPTION: It uses the following material models:
! DESCRIPTION:
! DESCRIPTION: 	- Mohr-Coulomb
! DESCRIPTION:
! DESCRIPTION:  - Hyperbolic Extended Drucker-Prager
!
! ELEMENT TYPE: SOLID45,SOLID185,SURF154
! MODULES:
! UNITS: SI
! KEYWORD1: Materials
! KEYWORD2: Geotechnics
! KEYWORD3: Material Models
!
!***********************************************************************
  FINISH
  ~CFCLEAR,,1
  NomFile='CFVR8003'
  /TITLE, %NomFile%, CivilFEM's Mohr-Coulomb

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

  *AFUN,DEG

  HYP   = 100.   ! hyperbolic parameter
  ECC   = 1.     ! eccentricity parameter
  LENG  = 1.     ! length of the cube edge
  NDIV  = 4      ! number of division per side
  NTINC = 50     ! number of substeps per step
  FS    = 1.5E+4 ! surface force
  INODE = 1      ! node for results checking

! ----------------------------------------------------------------------
! Model definition
! ----------------------------------------------------------------------
  LOCAL,11,0

  ! Materials

  ! Mohr-Coulomb
  ~CFMP,1,LIB,SOIL,,CL
  ~CFMP,1,SOIL,KPLA,,2
  ~CFMP,1,SOIL,KMCSP,,1
  ~CFMP,1,SOIL,HYP,,HYP
  ~CFMP,1,SOIL,ECC,,ECC
  ~CFMP,1,SOIL,IFLOW,,0
  ~CFGET,PHI,MATERIAL,1,SOIL,PHIMCeff
  ~CFGET,COH,MATERIAL,1,SOIL,CMCeff

  ! Hyperbolic Extended Drucker-Prager
  ~CFMP,2,LIB,SOIL,,CL
  TB,EDP,2,,,HYFUN
  TMP  = SIN(PHI)
  TMP1 = 6/(3 - TMP)
  C1 = TMP*TMP1
  C2 = C1*HYP
  C3 = COS(PHI)*TMP1*COH
  TBDATA,1,C1,C2,C3
  TB,EDP,2,,,HFPOT
  TBDATA,1,C1,C2

  ! Element types

  ET,1,SOLID45
  KEYOPT,1,1,1
  KEYOPT,1,2,0

  ET,2,SOLID185
  KEYOPT,2,2,0

  ET,3,SURF154
  KEYOPT,3,2,1
  KEYOPT,3,4,1

  ! Geometry

  K,1,0   ,   0,   0
  K,2,LENG,   0,   0
  K,3,LENG,LENG,   0
  K,4,   0,LENG,   0
  K,5,   0,   0,LENG
  K,6,LENG,   0,LENG
  K,7,LENG,LENG,LENG
  K,8,   0,LENG,LENG
  V,1,2,3,4,5,6,7,8
  LSEL,S,,,1,4,1
  LESIZE,ALL,,,NDIV

  ! Meshing

  LSEL,ALL
  TYPE,1
  EXTOPT,ESIZE,1
  VSWEEP,1,1,6
  TYPE,3
  ESYS,11
  AMESH,3

  ! Loads

  DA,1,UZ
  DA,6,UZ
  DA,5,UX
  DA,2,UY
  SFA,4,,PRES,4*FS
  SFA,2,,PRES,FS
  !ESEL,S,TYPE,,3
  !SFE,ALL,2,PRES,,FS

!-----------------------------------------------------------------------
! DATA CHECK
!-----------------------------------------------------------------------
! Data comparison number
  NCOMP    = 9
  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) = 'SEQV'
  LABEL(4) = 'EPELX'
  LABEL(5) = 'EPELY'
  LABEL(6) = 'EPELEQV'
  LABEL(7) = 'EPPLX'
  LABEL(8) = 'EPPLY'
  LABEL(9) = 'EPPLEQV'

! Warning and error tolerances
  TOLER(1,1) = 1.E-3 $ TOLER(1,2) = 1.E-3
  TOLER(2,1) = 1.E+3 $ TOLER(2,2) = 1.E+3
  TOLER(3,1) = 1.E+3 $ TOLER(3,2) = 1.E+3
  TOLER(4,1) = 1.E-5 $ TOLER(4,2) = 1.E-5
  TOLER(5,1) = 1.E-5 $ TOLER(5,2) = 1.E-5
  TOLER(6,1) = 1.E-5 $ TOLER(6,2) = 1.E-5
  TOLER(7,1) = 1.E-5 $ TOLER(7,2) = 1.E-5
  TOLER(8,1) = 1.E-5 $ TOLER(8,2) = 1.E-5
  TOLER(9,1) = 1.E-5 $ TOLER(9,2) = 1.E-5

! Correct values
!-----------------------------------------------------------------------
  /SOLU

  SOLCONTROL,OFF
  NLGEOM,OFF
  NROPT,FULL
  NEQIT,50
  NSUBST,NTINC

  OUTRES,ALL,NONE
  OUTRES,NSOL,LAST
  OUTRES,ESOL,LAST

  /PREP7
  ESLV
  EMODIF,ALL,MAT,2
  EMODIF,ALL,TYPE,2
  ESEL,ALL

  /SOLU
  SOLVE

  *GET,VALUE(1,1),NODE,INODE,S,X
  *GET,VALUE(2,1),NODE,INODE,S,Y
  *GET,VALUE(3,1),NODE,INODE,S,EQV
  *GET,VALUE(4,1),NODE,INODE,EPEL,X
  *GET,VALUE(5,1),NODE,INODE,EPEL,Y
  *GET,VALUE(6,1),NODE,INODE,EPEL,EQV
  *GET,VALUE(7,1),NODE,INODE,EPPL,X
  *GET,VALUE(8,1),NODE,INODE,EPPL,Y
  *GET,VALUE(9,1),NODE,INODE,EPPL,EQV

! Obtained values
!-----------------------------------------------------------------------
  /PREP7
  ESLV
  EMODIF,ALL,MAT,1
  EMODIF,ALL,TYPE,1
  ESEL,ALL

  /SOLU
  SOLVE

  *GET,VALUE(1,2),NODE,INODE,S,X
  *GET,VALUE(2,2),NODE,INODE,S,Y
  *GET,VALUE(3,2),NODE,INODE,S,EQV
  *GET,VALUE(4,2),NODE,INODE,EPEL,X
  *GET,VALUE(5,2),NODE,INODE,EPEL,Y
  *GET,VALUE(6,2),NODE,INODE,EPEL,EQV
  *GET,VALUE(7,2),NODE,INODE,EPPL,X
  *GET,VALUE(8,2),NODE,INODE,EPPL,Y
  *GET,VALUE(9,2),NODE,INODE,EPPL,EQV

!-----------------------------------------------------------------------
! Results comparison
!-----------------------------------------------------------------------
  COMPARA.MAC