!***********************************************************************
! TITLE: CFVR8001 - Mohr-Coulomb
! SUBTITLE: CivilFEM's Mohr-Coulomb equivalent to ANSYS's Linear Extended Drucker-Prager
!
! DESCRIPTION: Checking equivalence between the material models Mohr-Coulomb of
! DESCRIPTION: CivilFEM and Linear Extended Drucker-Prager of ANSYS
! DESCRIPTION:
! DESCRIPTION: Modelled acceleration is 9.81
! DESCRIPTION:
! DESCRIPTION: It uses the following material models:
! DESCRIPTION:
! DESCRIPTION: 	- Mohr-Coulomb
! DESCRIPTION:
! DESCRIPTION:  - Linear Extended Drucker-Prager
!
! ELEMENT TYPE: PLANE42,PLANE182,SURF153
! MODULES:
! UNITS: SI
! KEYWORD1: Materials
! KEYWORD2: Geotechnics
! KEYWORD3: Material Models
!
!***********************************************************************
  FINISH
  ~CFCLEAR,,1
  NomFile='CFVR8001'
  /TITLE, %NomFile%, CivilFEM's Mohr-Coulomb

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

  *AFUN,DEG

  HYP   = 0.  ! hyperbolic parameter
  ECC   = 1.  ! eccentricity parameter
  NSTAG = 1   ! number of phases
  H     = 40. ! slop heigth
  LENG  = 60. ! base length
  BET   = 35. ! slope inclination angle (0º < BET < 90º)
  NDIV  = 20  ! number of slope vertical divisions
  NTINC = 20  ! number of substeps per step
  INODE = 33  ! 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

  ! Linear Extended Drucker-Prager
  ~CFMP,2,LIB,SOIL,,CL
  TB,EDP,2,,,LYFUN
  TMP  = SIN(PHI)
  TMP1 = 6/(3 - TMP)
  C1 = TMP*TMP1
  C2 = COS(PHI)*TMP1*COH
  TBDATA,1,C1,C2
  TB,EDP,2,,,LFPOT
  TBDATA,1,C1

  ! Element types

  ET,1,PLANE42
  KEYOPT,1,2,1
  KEYOPT,1,3,2

  ET,2,PLANE182
  KEYOPT,2,1,0
  KEYOPT,2,3,2

  ET,3,SURF153
  KEYOPT,3,2,1
  KEYOPT,3,3,2
  KEYOPT,3,4,1

  ! Geometry

  DELY = H/NSTAG
  DELX = DELY/TAN(BET)
  K,1,0,0
  J = NSTAG + 2
  K,J,LENG,0
  *DO,I,1,NSTAG
    I = I + 1
    J = J + 1
    TMP = KY(I-1) + DELY
    K,I,0,TMP
    K,J,KX(J-1) - DELX,TMP
  *ENDDO
  TMP = NSTAG + 1
  *DO,I,1,NSTAG
    J = I + TMP
    A,I,J,J+1,I+1
  *ENDDO

  ! Meshing

  MSHKEY,0
  MSHAPE,1
  SNB  = SIN(BET)
  TMP  = NDIV/NSTAG
  TMP1 = TMP/SNB
  J = 4
  *DO,I,1,NSTAG
    LESIZE,J,,,TMP
    LESIZE,J-2,,,TMP1
    J = J + 3
  *ENDDO
  LESIZE,1,,,NDIV*LENG/H
  J = 3
  *DO,I,1,NSTAG
    LESIZE,J,,,NDIV*(LENG - I*DELX)/H
    J = J + 3
  *ENDDO

  ! Components

  *DO,I,1,NSTAG       ! Elements
    ASEL,S,AREA,,I
    AMESH,I
    ESLA,S
    CM,STAGE_%I%,ELEM
  *ENDDO

  LSEL,S,LINE,,1      ! Restrained nodes
  NSLL,S,1
  CM,NSET_ALL,NODE
  NSEL,NONE
  J = 4
  *DO,I,1,NSTAG
    LSEL,S,LINE,,J
    NSLL,A,1
    J = J + 3
  *ENDDO
  CM,NSET_UX,NODE

  ! Loads

  NSEL,ALL
  D,ALL,ALL

  ACEL,,9.81

!-----------------------------------------------------------------------
! 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.0E+4 $ TOLER(1,2) = 1.0E+4
  TOLER(2,1) = 1.0E+4 $ TOLER(2,2) = 1.0E+4
  TOLER(3,1) = 1.0E+4 $ TOLER(3,2) = 1.0E+4
  TOLER(4,1) = 1.3E-4 $ TOLER(4,2) = 1.3E-4
  TOLER(5,1) = 1.3E-4 $ TOLER(5,2) = 1.3E-4
  TOLER(6,1) = 1.0E-3 $ TOLER(6,2) = 1.0E-3
  TOLER(7,1) = 0.5E-4 $ TOLER(7,2) = 0.5E-4
  TOLER(8,1) = 0.5E-4 $ TOLER(8,2) = 0.5E-4
  TOLER(9,1) = 0.5E-4 $ TOLER(9,2) = 0.5E-4

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

  SOLCONTROL,OFF
  NLGEOM,ON
  SSTIF,ON
  NROPT,FULL
  NEQIT,50
  LNSRCH,ON
  NSUBST,NTINC
  ESTIF,1.E-08

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

  NSEL,ALL
  ESEL,ALL
  EKILL,ALL

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

  /SOLU
  *DO,I,1,NSTAG
    TIME,I

    CMSEL,S,STAGE_%I%
    EALIVE,ALL
    NSLE,S
    DDELE,ALL,ALL
    CMSEL,R,NSET_UX
    D,ALL,UX
    *IF,I,EQ,1,THEN
      NSLE,S
      CMSEL,R,NSET_ALL
      D,ALL,ALL
    *ENDIF

    ESEL,ALL
    NSEL,ALL
    SOLVE
  *ENDDO

  *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
  ESLA
  EMODIF,ALL,MAT,1
  EMODIF,ALL,TYPE,1

  /SOLU
  *DO,I,1,NSTAG
    TIME,I

    CMSEL,S,STAGE_%I%
    EALIVE,ALL
    NSLE,S
    DDELE,ALL,ALL
    CMSEL,R,NSET_UX
    D,ALL,UX
    *IF,I,EQ,1,THEN
      NSLE,S
      CMSEL,R,NSET_ALL
      D,ALL,ALL
    *ENDIF

    ESEL,ALL
    NSEL,ALL
    SOLVE
  *ENDDO

  *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