!*************************************************************************************!
! TITLE: CFVR1705 - Eurocode 8-1994: Seismic Analysis (Elastic spectrum) Y-axis
! SUBTITLE: Seismic Analysis on a 1 DOF beam according to Eurocode 8-1994
!
! DESCRIPTION: This example compares spectrum accelerations given by CivilFEM to accelerations
! DESCRIPTION: obtained when solving on a one degree of freedom beam and to the correct Eurocode spectum accelerations.
! DESCRIPTION:
! DESCRIPTION: The model is a single beam element with only one degree of freedom
! DESCRIPTION: (the rest of the DOF are constrained).
! DESCRIPTION: The beam has no mass. A punctual mass is applied at the free end of the beam.
! DESCRIPTION:
! DESCRIPTION: After solving, the mass acceleration is obtained and compared to the CivilFEM spectrum
! DESCRIPTION: acceleration and to the correct value. Therefore, this example checks the following results:
! DESCRIPTION:
! DESCRIPTION: - CivilFEM spectrum is defined correctly by comparing to the correct value.
! DESCRIPTION:
! DESCRIPTION: - CivilFEM spectrum is applied correctly and response acceleration is equal to spectrum acceleration.
!
! ELEMENT TYPE: BEAM4, MASS21
! MODULES:
! UNITS: User
! KEYWORD1: Seismic
! KEYWORD2: Eurocode 8
!
!*************************************************************************************!
   FINISH
  ~CFCLEAR,,1
  NomFile='CFVR1705'
  /TITLE, %NomFile%, Eurocode 8-1994: Seismic Analysis (Elastic spectrum) Y-axis

! ----------------------------------------------------------------------
! Model definition
! ----------------------------------------------------------------------
! Setup: Code & Units
~UNITS,,LENG,CM
~UNITS,,TIME,S
~UNITS,,FORC,KP
~CODESEL,EC3-92,EC2-91,,,EC8-94
/PREP7
! Preprocessor
! ----------------------------------------------------------------------
! Materials: A-42
 ~CFMP,1,LIB,STEEL,EA,A42
! Modify density (Rho = 0)
 ~CFMP,1,USER
 ~CFMP,1,DatGen,RHO,,0
! Element Types
  ET,1,BEAM4   ! Type 1: 3D Beam
  ET,2,MASS21  ! Type 2: Mass
! Sections
 ~SSECLIB,1,1,1,1   ! IPE 80
 ~BMSHPRO,1,BEAM,1,1,,,4,,0,,

! Nodes
  L = 10              ! L    : Bar length
  N, 1
  N,10,L

! Elements
  TYPE,1
  MAT ,1
  REAL,1
  EN,1,1,10
  EPLOT
! Parameters
  *GET,EXX,EX,1                         ! Ex   : Elastic Modulus
  ~CFGET,BarAr,SECTION,1,MECHPROP,A,,2  ! Area : Bar area
  ~CFGET,IZZ,SECTION,1,MECHPROP,IZZ,,2  ! Izz  : Moment of inertia
  ~CFGET,IYY,SECTION,1,MECHPROP,IYY,,2  ! Iyy  : Moment of inertia

! Spectrum definition
  ~DEFSPEC,ALL,1.2753,A,ELASTIC,1,1,5,

  pi=3.141592654
  Direction= 2 ! Direction ( 1: X, 2: Y, 3: Z),1

! ----------------------------------------------------------------------
! DATA CHECK
! ----------------------------------------------------------------------
! Data comparison number
  NComp = 40
  NComp_ch = 0

! Matrix dim.
  *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

*DO,II,1,20
   /PREP7
   ~CFGET,Per,SEISM,,SPEC,TH,,II,2   ! Period
   *IF,Direction,EQ,1,THEN
    ! Labels
     LABEL(II,1) = 'Sa(%II%)x'
     LABEL(20+II,1) = 'Sa(%II%)x'
     K=EXX*BarAr/L
     Mass=K*(Per/(2*pi))**2  ! Mass applied at end
     R,2,Mass
   *ELSEIF,Direction,EQ,2,THEN
    ! Labels
     LABEL(II,1) = 'Sa(%II%)y'
     LABEL(20+II,1) = 'Sa(%II%)y'
     K=12*EXX*IZZ/L**3      ! Bending stiffness around Z axis
     Mass=K*(Per/(2*pi))**2  ! Mass applied at end
     R,2,0,Mass
   *ELSEIF,Direction,EQ,3,THEN
    ! Labels
     LABEL(II,1) = 'Sa(%II%)z'
     LABEL(20+II,1) = 'Sa(%II%)z'
     K=12*EXX*IYY/L**3      ! Bending stiffness around Y axis
     Mass=K*(Per/(2*pi))**2  ! Mass applied at end
     R,2,0,0,Mass
   *ENDIF
   TYPE,2
   MAT, 1
   REAL,2
   EN,2,10
   FINISH
   /DELETE,file,rcv
   /DELETE,file,rst
   /SOLU

 ! Solution
 ! ---------------------------------------------------------------------
 ! Displacements
   D,1,,,,,,ALL,,,,,
   *IF,Direction,EQ,1,THEN
     D,10,,,,,,UY,UZ,ROTX,ROTY,ROTZ,
   *ELSEIF,Direction,EQ,2,THEN
     D,10,,,,,,UX,UZ,ROTX,ROTY,ROTZ,
   *ELSEIF,Direction,EQ,3,THEN
     D,10,,,,,,UX,UY,ROTX,ROTY,ROTZ,
   *ENDIF
   ~MODLSOL,1
   *IF,Direction,EQ,1,THEN
     ~CMBMOD,HORIZONT,NONE,NONE,0,0.0001,
   *ELSEIF,Direction,EQ,2,THEN
     ~CMBMOD,NONE,HORIZONT,NONE,0,0.0001,
   *ELSEIF,Direction,EQ,3,THEN
     ~CMBMOD,NONE,NONE,HORIZONT,0,0.0001,
   *ENDIF
 ! Query results
   *IF,Direction,EQ,1,THEN
     *GET,Displ,NODE,10,U,X
   *ELSEIF,Direction,EQ,2,THEN
     *GET,Displ,NODE,10,U,Y
   *ELSEIF,Direction,EQ,3,THEN
     *GET,Displ,NODE,10,U,Z
   *ENDIF
 ! Correct values
   VALUE(20-II+1,1)=Displ*K/Mass ! Acceleration
 ! Obtained values
   ~CFGET,VALUE(20-II+1,2),SEISM,,SPEC,SDH,,II,2	 ! Y components of elastic horizontal spectrum
   ~CFGET,VALUE(40-II+1,2),SEISM,,SPEC,SDH,,II,2	 ! Y components of elastic horizontal spectrum
   ! CivilFEM gives the spectrum Normalized by gravity.
   ! To compare it must be changed to acceleration user units.
   VALUE(20-II+1,2) = VALUE(20-II+1,2)*981
   VALUE(40-II+1,2) = VALUE(40-II+1,2)*981
   /PREP7
   EDELE,2
*ENDDO

! Correct values
! ----------------------------------------------------------------------
	  VALUE(40,1) = 0.3825E-01
	  VALUE(39,1) = 0.4250
	  VALUE(38,1) = 0.7286
	  VALUE(37,1) = 1.133
	  VALUE(36,1) = 1.569
	  VALUE(35,1) = 1.943
	  VALUE(34,1) = 2.550
	  VALUE(33,1) = 2.833
	  VALUE(32,1) = 3.188
	  VALUE(31,1) = 3.188
	  VALUE(30,1) = 3.188
	  VALUE(29,1) = 3.188
	  VALUE(28,1) = 2.951
	  VALUE(27,1) = 2.714
	  VALUE(26,1) = 2.478
	  VALUE(25,1) = 2.241
	  VALUE(24,1) = 2.004
	  VALUE(23,1) = 1.767
	  VALUE(22,1) = 1.531
	  VALUE(21,1) = 1.294

! Warning and error tolerances
! ----------------------------------------------------------------------
  *DO,II,1,NComp
   TOLER(II, 1)= 1E-03 $ TOLER(II, 2)= 1E-03
  *ENDDO

! ----------------------------------------------------------------------
! Results Comparison
! ----------------------------------------------------------------------
  COMPARA.MAC