!*************************************************************************************!
! TITLE: CFVR5241 - Isostatic and hiperestatic checking on isostatic beam
! SUBTITLE: rectangular cross section.
!
! DESCRIPTION: The model is a prestressed beam with a length of 15 meters.
! DESCRIPTION:
! DESCRIPTION: The beam is made of C16/20 concrete and the tendons are made of Gr250 steel.
! DESCRIPTION:
! DESCRIPTION: The section is rectangular 1 x 0.5, with reinforcement of S400 steel.
! DESCRIPTION:
! DESCRIPTION: Load states:
! DESCRIPTION:
! DESCRIPTION: - LS 1: Self weight
! DESCRIPTION:
! DESCRIPTION: - LS 2: Self weight + Prestressing loads
! DESCRIPTION:
! DESCRIPTION: - LS 3: Transverse weight
! DESCRIPTION:
! DESCRIPTION: - LS 4: Prestressing loads
! DESCRIPTION:
! DESCRIPTION: The aim of the example is to verify the calculation of prestressing losses.
! DESCRIPTION:
! DESCRIPTION: This example is based on "Prestressed concrete, a fundamental approach",
! DESCRIPTION: 4<sup>th</sup> edition, page 16.
!
! ELEMENT TYPE: BEAM44
! MODULES: PRSC
! UNITS: SI
! KEYWORD1: Beam
! KEYWORD2: Tendon
! KEYWORD3: Cross Section
! KEYWORD4: Prestressed Concrete
! KEYWORD5: Prestressing Losses
!
!*************************************************************************************!
   FINISH
   ~CFCLEAR,,1
   ~CFACTIV,PRSC,Y
   NomFile='CFVR5241'
 /TITLE,%NomFile%, I&H CHECKING, ISOST BEAM, RECT SECTION

! -------------------------------------------------------------------------------------
! Initial data
! -------------------------------------------------------------------------------------
! CivilFEM SETUP
  ~UNITS,SI
  ~CODESEL,,EC2-91,EC2-91
  ~CFACTIV,PRSC,Y

! Preprocessor
  /PREP7

! Material Definition
  ~CFMP,1,LIB,CONCRETE,EC2,C16/20
  ~CFMP,2,LIB,REINF,CEB,S400
  ~CFMP,3,LIB,PREST,ASTMA416,Gr250

 ! Parameters to ignore prestressing losses.
  ~CFMP,3,Prest,MU,,0
  ~CFMP,3,Prest,K ,,0
  ~CFMP,3,Prest,A ,,0

! Element Type Definition
  ET,1,Beam44

! Cross-section definition
  ~CSECDMS,1,REC,1,1 ,0.5,0 ,0 ,0 ,0 ,0 ,0 ,1 ,2
! Beam & Shell property
  ~BMSHPRO,1,BEAM,1,1,,,44,1,1,,Beam 1

! Model Construction
  k,1
  k,2,15
  l,1,2
  dk,1,ux,0,,,uy,uz,rotx
  dk,2,uy,0,,,uz
  h_elems = 4
  medcut = h_elems+1
  endcut   = h_elems*2+1
  lesize,all,,,h_elems*2
  lmesh,all

! Support beam construction
  area = 0.01
  c_diameter = 0.15
  ~SBBMDEF

  ~PCEPDEF,1,       1,     0,0
  ~PCEPDEF,2,medcut,-0.400,0
  ~PCEPDEF,3,  endcut,     0,0

  ~PCPPDEF,1,       1,     0,0
  ~PCPPDEF,3,  endcut,     0,0

  Papp = 862207
  ~PCTNDEF,1,3,area,c_diameter,Papp  ,Papp  ,0

  ~PCTNMDF,1,EADD,1.
  ~PCTNMDF,1,EADD,2.
  ~PCTNMDF,1,EADD,3.
  ~PCTNMDF,1,PADD,1.
  ~PCTNMDF,1,PADD,3.
  ~PCTNMDF,1,METHOD,1

  ~PCTNMDF,1,EPAROUT,1,0
  ~PCTNMDF,1,EPARIN,2,0

  ~PCLOSS,0,0,0,0.8,

/SOLU
  ACEL,,9.81
  SOLVE               ! LS 1: Self weight

  ~PCPL               ! Prestressing forces
  SOLVE               ! LS 2: SW + Prestressing loads

  ~PCDEL              ! Prestressing forces delete
  ACEL,,,9.81         ! "Transverse" weight
  SOLVE               ! LS 3: SWT

  ACEL,
  ~PCPL               ! Prestressing forces
  SOLVE               ! LS 4: Prestressing loads

/POST1

!--------------------------------------------------------------------------------------
! DATA CHECK
!--------------------------------------------------------------------------------------
! Data comparison number

  NComp = 24
  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
  *DIM,NODOS,,Ncomp,1
  *DIM,COORDX,,Ncomp,1

! Labels, Correct data, obtained data
!--------------------------------------------------------------------------------------
  ~CFSET,,1
  LABEL( 1) = 'KPRES1' $  VALUE( 1,1) =  0.0           $  ~CFGET,VALUE( 1,2),PRSCONC,,KPRES
  ~CFSET,,2
  LABEL( 2) = 'KPRES2' $  VALUE( 2,1) =  1             $  ~CFGET,VALUE( 2,2),PRSCONC,,KPRES
                                   ~CFGET,aux1,ELEMENT,h_elems,FORCE,HMZ,J,
                                   ~CFGET,aux2,ELEMENT,h_elems,FORCE,IMZ,J,
  LABEL( 3) = 'HMZ_IMZ'$  VALUE( 3,1) = -1             $  VALUE( 3,2)=aux1/aux2
  LABEL( 4) = 'FY'     $  VALUE( 4,1) =  0.0           $  ~CFGET,VALUE( 4,2),ELEMENT,h_elems,FORCE, FY,J,
  LABEL( 5) = 'IFY'    $  VALUE( 5,1) =  0.0           $  ~CFGET,VALUE( 5,2),ELEMENT,h_elems,FORCE,IFY,J,

  ~CFSET,,3
  LABEL( 6) = 'KPRES3' $  VALUE( 6,1) =  0.0           $  ~CFGET,VALUE( 6,2),PRSCONC,,KPRES
  LABEL( 7) = 'MY'     $  VALUE( 7,1) =  -344882.8125  $  ~CFGET,VALUE( 7,2),ELEMENT,h_elems,FORCE, MY,J,
  LABEL( 8) = 'IMY'    $  VALUE( 8,1) =  0.0           $  ~CFGET,VALUE( 8,2),ELEMENT,h_elems,FORCE,IMY,J,

  ~CFSET,,4
  LABEL( 9) = 'KPRES4' $  VALUE( 9,1) =  1             $  ~CFGET,VALUE( 9,2),PRSCONC,,KPRES
                                   ~CFGET,aux1,ELEMENT,h_elems,FORCE,HMZ,J,
                                   ~CFGET,aux2,ELEMENT,h_elems,FORCE,IMZ,J,
  LABEL(10) = 'HMZ_IMZ'$  VALUE(10,1) =  0.0           $  VALUE(10,2)=ABS(aux1/aux2)
                                   ~CFGET,aux1,ELEMENT,h_elems,FORCE,HFX,J,
                                   ~CFGET,aux2,ELEMENT,h_elems,FORCE,IFX,J,
  LABEL(11) = 'HFX_IFX'$  VALUE(11,1) =  0.0           $  VALUE(11,2)=ABS(aux1/aux2)
                                   ~CFGET,aux1,ELEMENT,h_elems,FORCE,HFY,J,
                                   ~CFGET,aux2,ELEMENT,h_elems,FORCE,IFY,J,
  LABEL(12) = 'HFY'     $ VALUE(12,1) =  0.0           $  VALUE(12,2)=ABS(aux1)
  LABEL(13) = 'RFX1'    $ VALUE(13,1) =  0.0           $  *GET,VALUE(13,2),NODE,1,RF,FX
  LABEL(14) = 'RFY1'    $ VALUE(14,1) =  0.0           $  *GET,VALUE(14,2),NODE,1,RF,FY
  LABEL(15) = 'RFZ1'    $ VALUE(15,1) =  0.0           $  *GET,VALUE(15,2),NODE,1,RF,FZ
  LABEL(16) = 'RFY2'    $ VALUE(16,1) =  0.0           $  *GET,VALUE(16,2),NODE,2,RF,FY
  LABEL(17) = 'RFZ2'    $ VALUE(17,1) =  0.0           $  *GET,VALUE(17,2),NODE,2,RF,FZ

  a_parab = 0.0071111111
  LABEL(18) = 'YP'
  VALUE(18,1) = a_parab*(-7.5+7.5/h_elems)**2-0.4
  ~CFGET,VALUE(18,2),PRSCONC,,TENDON,COORD,Y,1,2
  LABEL(19) = 'ZP'
  VALUE(19,1) = 0.0
  ~CFGET,VALUE(19,2),PRSCONC,,TENDON,COORD,Z,1,2
  LABEL(20) = 'SXYP'
  VALUE(20,1) = 2*a_parab*(-7.5+7.5/h_elems)
  ~CFGET,VALUE(20,2),PRSCONC,,TENDON,COORD,SXY,1,2

  ~CFSET,,2
  LABEL(21) = 'UYMAX'
  VALUE(21,1) =  0.0
  *GET,VALUE(21,2),NODE, 1,U,Y
  LABEL(22) = 'UYMIN'
  VALUE(22,1) =  -0.498E-5
  *GET,VALUE(22,2),NODE, 6,U,Y
  LABEL(23) = 'UZMAX'
  VALUE(23,1) =  0.0
  *GET,VALUE(23,2),NODE, 1,U,Z
  LABEL(24) = 'UZMIN'
  VALUE(24,1) =  0.0
  *GET,VALUE(24,2),NODE, 6,U,Z

! Warning and error tolerances
  TOLER(  1,1) = 1E-02   $  TOLER(  1,2) = 1E-02 
  TOLER(  2,1) = 1E-02   $  TOLER(  2,2) = 1E-02 
  TOLER(  3,1) = 1E-02   $  TOLER(  3,2) = 1E-02   
  TOLER(  4,1) = 1E-02   $  TOLER(  4,2) = 1E-02 
  TOLER(  5,1) = 1E-02   $  TOLER(  5,2) = 1E-02 
  TOLER(  6,1) = 1E-02   $  TOLER(  6,2) = 1E-02 
  TOLER(  7,1) = 1E-00   $  TOLER(  7,2) = 1E-00 
  TOLER(  8,1) = 1E-02   $  TOLER(  8,2) = 1E-02 
  TOLER(  9,1) = 1E-02   $  TOLER(  9,2) = 1E-02 
  TOLER( 10,1) = 1E-02   $  TOLER( 10,2) = 1E-02 
  TOLER( 11,1) = 1E-02   $  TOLER( 11,2) = 1E-02 
  TOLER( 12,1) = 1E-02   $  TOLER( 12,2) = 1E-02 
  TOLER( 13,1) = 1E-02   $  TOLER( 13,2) = 1E-02
  TOLER( 14,1) = 1E-02   $  TOLER( 14,2) = 1E-02 
  TOLER( 15,1) = 1E-02   $  TOLER( 15,2) = 1E-02
  TOLER( 16,1) = 1E-02   $  TOLER( 16,2) = 1E-02 
  TOLER( 17,1) = 1E-02   $  TOLER( 17,2) = 1E-02 
  TOLER( 18,1) = 1E-02   $  TOLER( 18,2) = 1E-02 
  TOLER( 19,1) = 1E-02   $  TOLER( 19,2) = 1E-02 
  TOLER( 20,1) = 1E-02   $  TOLER( 20,2) = 1E-02 
  TOLER( 21,1) = 1E-02   $  TOLER( 21,2) = 1E-02 
  TOLER( 22,1) = 1E-08   $  TOLER( 22,2) = 1E-08 

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