/POST1
SET,LAST
*DEL,ALL

! ================== Safety Factors and Material Selection ====================
SF_FF = 2.223900 !Safety factor for fibre failure
SF_IFF = 0.832922 !Safety factor for inter-fibre failure
MatNo=0 !Material number to apply Puck to (set to 0 for worst index in ply stack
PlyAngle=0.000000 !Ply angle if material number not equal to 0
FailMode=3 !Ply angle if material number not equal to 0
Case=2 !Load case number
Factor=1 !Factor applied to loads

! ================ Script Starts ===========================
!Create material property arrays
*DIM,XTE,ARRAY,12,1,1, , ,
*DIM,XCE,ARRAY,12,1,1, , ,
*DIM,SS,ARRAY,12,1,1, , ,
*DIM,YT,ARRAY,12,1,1, , ,
*DIM,YC,ARRAY,12,1,1, , ,
*DIM,PUCKPLUS,ARRAY,12,1,1, , ,
*DIM,PUCKMINUS,ARRAY,12,1,1, , ,
*DIM,PUCK,ARRAY,12,1,1, , ,
! Puck properties for Gelcoat
*SET,XTE(1,1,1) , 0.000000/SF_FF
*SET,XCE(1,1,1) , 0.000000/SF_FF
*SET,SS(1,1,1) , 0.000000/SF_IFF
*SET,YT(1,1,1) , 0.000000/SF_IFF
*SET,YC(1,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(1,1,1) , 0.000000
*SET,PUCKMINUS(1,1,1) , 0.000000
*SET,PUCK(1,1,1) , 0.000000
! Puck properties for steel
*SET,XTE(2,1,1) , 0.000000/SF_FF
*SET,XCE(2,1,1) , 0.000000/SF_FF
*SET,SS(2,1,1) , 0.000000/SF_IFF
*SET,YT(2,1,1) , 0.000000/SF_IFF
*SET,YC(2,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(2,1,1) , 0.000000
*SET,PUCKMINUS(2,1,1) , 0.000000
*SET,PUCK(2,1,1) , 0.000000
! Puck properties for steel_drive
*SET,XTE(3,1,1) , 0.000000/SF_FF
*SET,XCE(3,1,1) , 0.000000/SF_FF
*SET,SS(3,1,1) , 0.000000/SF_IFF
*SET,YT(3,1,1) , 0.000000/SF_IFF
*SET,YC(3,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(3,1,1) , 0.000000
*SET,PUCKMINUS(3,1,1) , 0.000000
*SET,PUCK(3,1,1) , 0.000000
! Puck properties for cast_iron
*SET,XTE(4,1,1) , 0.000000/SF_FF
*SET,XCE(4,1,1) , 0.000000/SF_FF
*SET,SS(4,1,1) , 0.000000/SF_IFF
*SET,YT(4,1,1) , 0.000000/SF_IFF
*SET,YC(4,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(4,1,1) , 0.000000
*SET,PUCKMINUS(4,1,1) , 0.000000
*SET,PUCK(4,1,1) , 0.000000
! Puck properties for Glass_UD
*SET,XTE(5,1,1) , 0.017966/SF_FF
*SET,XCE(5,1,1) , -0.014988/SF_FF
*SET,SS(5,1,1) , 90300000.000000/SF_IFF
*SET,YT(5,1,1) , 48600000.000000/SF_IFF
*SET,YC(5,1,1) , 157000000.000000/SF_IFF
*SET,PUCKPLUS(5,1,1) , 0.300000
*SET,PUCKMINUS(5,1,1) , 0.250000
*SET,PUCK(5,1,1) , 1.000000
! Puck properties for Carbon_UD
*SET,XTE(6,1,1) , 0.015127/SF_FF
*SET,XCE(6,1,1) , -0.007481/SF_FF
*SET,SS(6,1,1) , 55000000.000000/SF_IFF
*SET,YT(6,1,1) , 46600000.000000/SF_IFF
*SET,YC(6,1,1) , 157700000.000000/SF_IFF
*SET,PUCKPLUS(6,1,1) , 0.350000
*SET,PUCKMINUS(6,1,1) , 0.300000
*SET,PUCK(6,1,1) , 1.000000
! Puck properties for Glass_biax
*SET,XTE(7,1,1) , 0.022388/SF_FF
*SET,XCE(7,1,1) , -0.013806/SF_FF
*SET,SS(7,1,1) , 144000000.000000/SF_IFF
*SET,YT(7,1,1) , 300000000.000000/SF_IFF
*SET,YC(7,1,1) , 185000000.000000/SF_IFF
*SET,PUCKPLUS(7,1,1) , 0.000000
*SET,PUCKMINUS(7,1,1) , 0.000000
*SET,PUCK(7,1,1) , 0.000000
! Puck properties for Glass_triax
*SET,XTE(8,1,1) , 0.024721/SF_FF
*SET,XCE(8,1,1) , -0.010623/SF_FF
*SET,SS(8,1,1) , 99300000.000000/SF_IFF
*SET,YT(8,1,1) , 127000000.000000/SF_IFF
*SET,YC(8,1,1) , 127000000.000000/SF_IFF
*SET,PUCKPLUS(8,1,1) , 0.300000
*SET,PUCKMINUS(8,1,1) , 0.250000
*SET,PUCK(8,1,1) , 1.000000
! Puck properties for medium_density_foam
*SET,XTE(9,1,1) , 0.000000/SF_FF
*SET,XCE(9,1,1) , 0.000000/SF_FF
*SET,SS(9,1,1) , 0.000000/SF_IFF
*SET,YT(9,1,1) , 0.000000/SF_IFF
*SET,YC(9,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(9,1,1) , 0.000000
*SET,PUCKMINUS(9,1,1) , 0.000000
*SET,PUCK(9,1,1) , 0.000000
! Puck properties for resin
*SET,XTE(10,1,1) , 0.000000/SF_FF
*SET,XCE(10,1,1) , 0.000000/SF_FF
*SET,SS(10,1,1) , 0.000000/SF_IFF
*SET,YT(10,1,1) , 0.000000/SF_IFF
*SET,YC(10,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(10,1,1) , 0.000000
*SET,PUCKMINUS(10,1,1) , 0.000000
*SET,PUCK(10,1,1) , 0.000000
! Puck properties for remove
*SET,XTE(11,1,1) , 0.000000/SF_FF
*SET,XCE(11,1,1) , 0.000000/SF_FF
*SET,SS(11,1,1) , 0.000000/SF_IFF
*SET,YT(11,1,1) , 0.000000/SF_IFF
*SET,YC(11,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(11,1,1) , 0.000000
*SET,PUCKMINUS(11,1,1) , 0.000000
*SET,PUCK(11,1,1) , 0.000000
! Puck properties for Glue
*SET,XTE(12,1,1) , 0.000000/SF_FF
*SET,XCE(12,1,1) , 0.000000/SF_FF
*SET,SS(12,1,1) , 0.000000/SF_IFF
*SET,YT(12,1,1) , 0.000000/SF_IFF
*SET,YC(12,1,1) , 0.000000/SF_IFF
*SET,PUCKPLUS(12,1,1) , 0.000000
*SET,PUCKMINUS(12,1,1) , 0.000000
*SET,PUCK(12,1,1) , 0.000000
!Get the highest number of layers
*GET,num_elems_,ELEM,0,COUNT !Get the number of elements
*GET,elem_,ELEM,0,NUM,MIN !Get label of the first element
MaxLay=0
*DO,i,1,num_elems_,1
	! Define some parameters
	*GET,Mat_,ELEM,elem_,ATTR,MAT
	*GET,Sect,ELEM,elem_,ATTR,SECN
    !Get the section details
	*GET, NLAYER, SHEL, Sect, Prop, NLAY !Get number of layers
	*IF,NLAYER,GT,MaxLay,THEN
		MaxLay=NLAYER
	*ENDIF
	*GET,elem_,ELEM,elem_,NXTH
*ENDDO
!Initialise the stress element tables
ETABLE,ERAS
!Create the output element tables
ETABLE,PuckFT,EPTH,X
ETABLE,PuckFTL,EPTH,Y
ETABLE,PuckFC,EPTH,Z
ETABLE,PuckFCL,EPTH,XY
ETABLE,PuckIFA,EPTH,YZ
ETABLE,PuckIFAL,EPTH,XZ
ETABLE,PuckIFB,EPTH,1
ETABLE,PuckIFBL,EPTH,2
ETABLE,PuckIFC,EPTH,3   
ETABLE,PuckIFCL,EPTH,INT
!Stress/strain variables
ETABLE,SX,S,X
ETABLE,SY,S,Y
ETABLE,SXY,S,XY
ETABLE,EX,EPEL,X
!Loop over the layers
Counter=0
/NOPR
*DO,i,1,MaxLay,1
	!Select this layer
	Layer,i
	RSYS,LSYS
	!Refill stress tables
	ETABLE,1,ERAS
	ETABLE,2,ERAS
	ETABLE,3,ERAS
	ETABLE,4,ERAS
	ETABLE,SX,S,X
	ETABLE,SY,S,Y
	ETABLE,SXY,S,XY
	ETABLE,EX,EPEL,X
	!Loop over the elements in this layer
	*GET,num_elems_,ELEM,0,COUNT !Get the number of elements
	*GET,elem_,ELEM,0,NUM,MIN !Get label of the first element
	*DO,j,1,num_elems_,1
		!Progess report
		*IF,Counter,EQ,1000,THEN
			Counter=0
 *MSG, INFO--,j,num_elems_,i,MaxLay
 Checking Element %6.0d of %6.0d in layer %4.0d of %4.0d 
		*ENDIF
		Counter=Counter+1
		! Define some parameters
		*GET,Mat_,ELEM,elem_,ATTR,MAT
		*GET,Sect,ELEM,elem_,ATTR,SECN
		!Get the section details
		*GET, NLAYER, SHEL, Sect, Prop, NLAY !Get number of layers
		!Check if this element has this layer
		*IF,i,LE,NLAYER,THEN
			!Check if this material needs to be analysed for Puck
			*GET, MATID, SHEL, Sect, LAYD, i, MAT
			*GET, AngleI, SHEL, Sect, LAYD, i, ANGL
			Analyse=0
			*IF,PUCK(MATID),EQ,1,THEN
				*IF,MatNo,EQ,0,THEN
					Analyse=1
				*ELSEIF,MATID,EQ,MATNO,THEN
					*IF,AngleI,EQ,PlyAngle,THEN
						Analyse=1
					*ENDIF
				*ENDIF
			*ENDIF
			!Which failure modes will be checked?
			*IF,FailMode,EQ,1,THEN
				FibreFailure=1
				InterFibreFailure=0
			*ELSEIF,FailMode,EQ,2,THEN
				FibreFailure=0
				InterFibreFailure=1
			*ELSEIF,FailMode,EQ,3,THEN
				FibreFailure=1
				InterFibreFailure=1
			*ENDIF
			*IF,ANALYSE,EQ,1,THEN
				!We analyse this element for Puck
				*GET,S1,ETAB,11,ELEM,elem_
				*GET,S2,ETAB,12,ELEM,elem_
				*GET,S12,ETAB,13,ELEM,elem_
				*GET,E1,ETAB,14,ELEM,elem_
				*GET,FT,ETAB,1,ELEM,elem_
				*GET,FC,ETAB,3,ELEM,elem_
				*GET,IFA,ETAB,5,ELEM,elem_
				*GET,IFB,ETAB,7,ELEM,elem_
				*GET,IFC,ETAB,9,ELEM,elem_
				*IF,FibreFailure,EQ,1,THEN
				!Check for fibre failure
				*IF,E1,GT,0,THEN
					!Tensile failure
					*IF,E1/XTE(MATID),GT,FT,THEN
						DETAB,elem_,PuckFT,E1/XTE(MATID)
						DETAB,elem_,PuckFTL,i
					*ENDIF
				*ENDIF
				*IF,E1,LE,0,THEN
					!Compressive failure
					*IF,E1/XCE(MATID),GT,FC,THEN
						DETAB,elem_,PuckFC,E1/XCE(MATID)
						DETAB,elem_,PuckFCL,i
					*ENDIF
				*ENDIF
				*ENDIF
				*IF,InterFibreFailure,EQ,1,THEN
				!Check for matrix failure A
				*IF,S2,GT,0,THEN
					FMA=((S12/SS(MATID))**2+(1-PUCKPLUS(MATID)*YT(MATID)/SS(MATID))**2*(S2/YT(MATID))**2)**0.5+PUCKPLUS(MATID)*S2/SS(MATID)
					*IF,FMA,GT,IFA,THEN
						DETAB,elem_,PuckIFA,FMA
						DETAB,elem_,PuckIFAL,i
					*ENDIF
				*ENDIF
				!Check for matrix failure B and C
				RA=SS(MATID)/(2*PUCKMINUS(MATID))*((1+2*PUCKMINUS(MATID)*YC(MATID)/SS(MATID))**0.5-1)
				PMNN=PUCKMINUS(MATID)*RA/SS(MATID)
				T21C=SS(MATID)*(1+2*PMNN)**0.5
				*IF,abs(S2/S12),LE,RA/ABS(T21C),THEN
					FMB=1/SS(MATID)*((S12**2+(PUCKMINUS(MATID)*S2)**2)**0.5+PUCKMINUS(MATID)*S2)
					*IF,FMB,GT,IFB,THEN
						DETAB,elem_,PuckIFB,FMB
						DETAB,elem_,PuckIFBL,i
					*ENDIF
				*ENDIF
				*IF,abs(S12/S2),LE,ABS(T21C)/RA,THEN
					FMC=((S12/(2*(1+PMNN)*SS(MATID)))**2+(S2/YC(MATID))**2)*(YC(MATID)/(-S2))
					*IF,FMC,GT,IFC,THEN
						DETAB,elem_,PuckIFC,FMC
						DETAB,elem_,PuckIFCL,i
					*ENDIF
				*ENDIF
				*ENDIF
			*ENDIF
		*ENDIF
		*GET,elem_,ELEM,elem_,NXTH
	*ENDDO
*ENDDO
SAVE, , , , ALL
! Write outputs 
*GET,num_elems_,ELEM,0,COUNT !Get the number of elements
*DIM,OutArray,ARRAY,num_elems_,10,1, , ,
*GET,elem_,ELEM,0,NUM,MIN !Get label of the first element
*DO,i,1,num_elems_,1
	*GET,FT,ETAB,1,ELEM,elem_
	*GET,FC,ETAB,3,ELEM,elem_
	*GET,IFA,ETAB,5,ELEM,elem_
	*GET,IFB,ETAB,7,ELEM,elem_
	*GET,IFC,ETAB,9,ELEM,elem_
	*GET,N1_,ELEM,elem_,NODE,1
	*GET,N2_,ELEM,elem_,NODE,2
	*GET,N3_,ELEM,elem_,NODE,3
	*GET,N4_,ELEM,elem_,NODE,4
	*SET,OutArray(i,1,1) , elem_
	*SET,OutArray(i,2,1) , N1_
	*SET,OutArray(i,3,1) , N2_
	*SET,OutArray(i,4,1) , N3_
	*SET,OutArray(i,5,1) , N4_
	*SET,OutArray(i,6,1) , FT
	*SET,OutArray(i,7,1) , FC
	*SET,OutArray(i,8,1) , IFA
	*SET,OutArray(i,9,1) , IFB
	*SET,OutArray(i,10,1) , IFC
	*GET,elem_,ELEM,elem_,NXTH
*ENDDO
*MWRITE, OutArray,FailureCriteria_Case-%Case%_Factor-%Factor%_Material-%MatNo%_PlyAngle-%PlyAngle%_Mode-%FailMode%,txt
(,F8.0,',',F8.0,',',F8.0,',',F8.0,',',F8.0,',',e16.8,',',e16.8,',',e16.8,',',e16.8,',',e16.8)
FINISH