!*******************************************************************************
! TITLE: CFVR3704 - Initial stresses on 3D Model
! SUBTITLE: Initial stresses on a 3D Model with irregular mesh.
!
! DESCRIPTION: A prism model is free meshed with tetrahedral-shaped elements.  
! DESCRIPTION: The model is rotated about a local coordinate system.
! DESCRIPTION: Initial stresses and gravity are applied.
! DESCRIPTION:
!
! ELEMENT TYPE: SOLID185
! MODULES: GETC
! UNITS: SI
! KEYWORD1: Geotechnics
!*******************************************************************************

FINISH
~CFCLEAR,,1
~CFACTIV,GETC,Y
~UNITS,SI
*AFUN,DEG

/DELETE,cfvr3704,ist
NomFile='CFVR3704'
/TITLE,%NomFile%,Initial stresses. 3D Model
 
! parameters
! ----------
RXY  = 30. ! Euler angles
RYZ  = -20.
RZX  = 55.
DEPT = 1. ! depth of the prism
NDIV = 30 ! number of divisions of the lines parallel to the axis Z 
NEL  = 4  ! number of elements where to check values
 
! element types
! -------------
/PREP7
ET,1,SOLID185

! materials
! ---------
~CFMP,1,LIB,SOIL,,GC
~CFGET,GAM,MATERIAL,1,DATGEN,GAM
~CFGET,K0,MATERIAL,1,SOIL,K0
 
! solid model and mesh
! --------------------
ESIZ = DEPT/NDIV
CLOCAL,11,0,,,,RXY,RYZ,RZX
TMP = 5*ESIZ
K,1,0.,0.,0.
K,2,TMP,0.,0.
K,3,TMP,TMP,0.
K,4,0.,TMP,0.
K,5,0.,0.,DEPT
K,6,TMP,0.,DEPT
K,7,TMP,TMP,DEPT
K,8,0.,TMP,DEPT

V,1,2,3,4,5,6,7,8
MSHAPE,1,3D
LESIZE,ALL,ESIZ
VMESH,ALL
NROTAT,ALL

! DATA CHECK
! ----------

NCOMP= NEL*6
NCOMP_CH = 0
*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

! rotation matrix
! ---------------
*DIM,ROT,ARRAY,3,3
CSXY      = COS(RXY)
SNXY      = SIN(RXY)
CSYZ      = COS(RYZ)
SNYZ      = SIN(RYZ)
CSZX      = COS(RZX)
SNZX      = SIN(RZX)
CSXY_SNYZ = CSXY*SNYZ
SNXY_SNYZ = SNXY*SNYZ
ROT(1,1)  =  CSXY*CSZX - SNXY_SNYZ*SNZX
ROT(2,1)  =  SNXY*CSZX + CSXY_SNYZ*SNZX
ROT(3,1)  = -CSYZ*SNZX
ROT(1,2)  = -SNXY*CSYZ
ROT(2,2)  =  CSXY*CSYZ
ROT(3,2)  =  SNYZ
ROT(1,3)  =  CSXY*SNZX + SNXY_SNYZ*CSZX
ROT(2,3)  =  SNXY*SNZX - CSXY_SNYZ*CSZX
ROT(3,3)  =  CSYZ*CSZX

! correct data 
! ------------
*DIM,IEL,,NEL ! elements where to check values
*DIM,SIG_L,,3,NEL
IEL(1) = 1,5,624,891

J = 1
*DO,I,1,NEL
  SIG_L(3,I) = -GAM*(DEPT - (CENTRX(IEL(I))*ROT(1,3) + CENTRY(IEL(I))*ROT(2,3) + CENTRZ(IEL(I))*ROT(3,3)))
  SIG_L(1,I) = K0*SIG_L(3,I)
  SIG_L(2,I) = SIG_L(1,I) 

  TMP  = SIG_L(1,I)*ROT(1,1)
  TMP1 = SIG_L(2,I)*ROT(1,2)
  TMP2 = SIG_L(3,I)*ROT(1,3)
  VALUE(J,1)   = ROT(1,1)*TMP + ROT(1,2)*TMP1 + ROT(1,3)*TMP2 ! SXX
  VALUE(J+3,1) = ROT(2,1)*TMP + ROT(2,2)*TMP1 + ROT(2,3)*TMP2 ! SXY
  VALUE(J+5,1) = ROT(3,1)*TMP + ROT(3,2)*TMP1 + ROT(3,3)*TMP2 ! SXZ
  TMP  = SIG_L(1,I)*ROT(2,1)
  TMP1 = SIG_L(2,I)*ROT(2,2)
  TMP2 = SIG_L(3,I)*ROT(2,3)
  VALUE(J+1,1) = ROT(2,1)*TMP + ROT(2,2)*TMP1 + ROT(2,3)*TMP2 ! SYY
  VALUE(J+4,1) = ROT(3,1)*TMP + ROT(3,2)*TMP1 + ROT(3,3)*TMP2 ! SYZ
  VALUE(J+2,1) = ROT(3,1)*SIG_L(1,I)*ROT(3,1) + ROT(3,2)*SIG_L(2,I)*ROT(3,2) + ROT(3,3)*SIG_L(3,I)*ROT(3,3) ! SZZ
  J = J + 6
*ENDDO

! obtained data 
! -------------
/SOLU

DA,1,UZ
DA,2,UY
DA,5,UX
SFGRAD,PRES,11,Z,DEPT,-GAM*K0
SFA,3,,PRES,0.
SFA,4,,PRES,0.
~TIS,cfvr3704,ist,Z
ACEL,9.8*ROT(1,3),9.8*ROT(2,3),9.8*ROT(3,3)
INISTATE,READ,cfvr3704,ist

ERESX,NO
NROPT,FULL
SOLVE

/POST1
SET,LAST

ETABLE,,S,X
ETABLE,,S,Y
ETABLE,,S,Z
ETABLE,,S,XY
ETABLE,,S,YZ
ETABLE,,S,XZ
J = 1
*DO,I,1,NEL
  *GET,VALUE(J  ,2),ETAB,1,ELEM,IEL(I)
  *GET,VALUE(J+1,2),ETAB,2,ELEM,IEL(I)
  *GET,VALUE(J+2,2),ETAB,3,ELEM,IEL(I)
  *GET,VALUE(J+3,2),ETAB,4,ELEM,IEL(I)
  *GET,VALUE(J+4,2),ETAB,5,ELEM,IEL(I)
  *GET,VALUE(J+5,2),ETAB,6,ELEM,IEL(I)
  J = J + 6
*ENDDO

! results comparison
! ----------------------------------------------------------------------
EREL = 1.0
*DO,I,1,NCOMP
  TMP = EREL*ABS(VALUE(I,2) - VALUE(I,1))
  TOLER(I,1) = TMP 
  TOLER(I,2) = TMP
*ENDDO

COMPARA.MAC