!MSc Thesis Code
!Author: Amin Deldari Alamdari
!Topic: Optimization of rib reinforced I-beams
!Department of Mechanical Engineering, Bogazici University

!-------------------------------------------------------------
!Units: m, kg, m/s^2, kg/m^3, N, Pa
!-------------------------------------------------------------!------------------------------------------
FINISH
/CLEAR,NOSTART

!*******************************************************************************************************
!								Initial OF THE PROGRAM
!*******************************************************************************************************
*CREATE,Initial

KEYW,PR_SGVOF,0 			!NO WARNING MESSAGE
/NERR,0,100000,,0,0 		!LIMITING THE NUMBER OF WARNING AND ERROR MESSAGES DISPLAYED
!/FTYPE,ALL,INT 
/PREP7
/NOPR   
OUTRES,ALL,NONE			!CONTROLLING THE SOLUTION DATA WRITTEN TO THE DATABASE
OUTRES,MISC
OUTRES,ESOL
OUTPR,ALL,LAST 
!/FILNAME,buckling,0
KEYW,PR_SET,1   
KEYW,PR_STRUC,1 
KEYW,PR_THERM,0 
KEYW,PR_FLUID,0 
KEYW,PR_ELMAG,0 
KEYW,MAGNOD,0   
KEYW,MAGEDG,0   
KEYW,MAGHFE,0   
KEYW,MAGELC,0   
KEYW,PR_MULTI,0 
!*
/VIEW,1,1,1,1   
/ANG,1  
/REP,FAST  
!*
*AFUN,DEG

!*******************************************************************************************************
!										ELEMENT TYPE
!*******************************************************************************************************


ET,1,SHELL181				!ELEMENT TYPE

KEYOPT,1,1,0 
KEYOPT,1,3,0 
KEYOPT,1,5,0 
KEYOPT,1,8,0 
KEYOPT,1,9,0 

!*******************************************************************************************************
!										MATERIAL PROPERTIES
!*******************************************************************************************************

ThicknessPly=0.000125				!PLY THICKNESS
LaminatePlyNumberTop=8				!PLY NUMBER OF THE TOP FLANGE LAMINATES
LaminatePlyNumberWeb=4				!PLY NUMBER OF THE WEB LAMINATES
LaminatePlyNumberBottom=8			!PLY NUMBER OF THE BOTTOM FLANGE LAMINATES

ElasticModulus1=165e9				!LONGITUDINAL ELASTIC MODULUS (E1) (FIBER DIRECTION)
ElasticModulus2=9e9					!TRANSVERSE ELASTIC MODULUS (E2)
ElasticModulus3=9e9					!THROUGH-THICKNESS ELASTIC MODULUS (E3)

ShearModulus12=5.6e9				!IN-PLANE SHEAR MODULUS (G12)
ShearModulus13=5.6e9				!TRANSVERSE SHEAR MODULUS (G13)
ShearModulus23=2.8e9				!THROUGH-THICKNESS SHEAR MODULUS (G23)

PoissonsRatio12=0.34				!MAJOR POISSON'S RATIO (v12)
PoissonsRatio13=0.34				!MAJOR TRANSVERSE POISSON'S RATIO (v13)
PoissonsRatio23=0.5					!THROUGH-THICKNESS POISSON'S RATIO (v23)

TensileStrength1=2560e6				!LONGITUDINAL TENSILE STRENGTH (XT) (FIBER DIRECTION)
CompressiveStrength1=-1590e6		!LONGITUDINAL COMPRESSIVE STRENGTH (XC)
TensileStrength2=73e6				!TRANSVERSE TENSILE STRENGTH (YT)
CompressiveStrength2=-185e6			!TRANSVERSE COMPRESSIVE STRENGTH (YC)
TensileStrength3=63e6				!THROUGH-THICKNESS TENSILE STRENGTH (ZT)
CompressiveStrength3=-185e6			!THROUGH-THICKNESS COMPRESSIVE STRENGTH (ZC)

ShearStrength12=90e6				!IN-PLANE SHEAR STRENGTH (S12)
ShearStrength13=90e6				!TRANSVERSE SHEAR STRENGTH (S13)
ShearStrength23=57e6				!THROUGH-THICKNESS SHEAR STRENGTH (S23)

TensileStrain1=0.01551				!LONGITUDINAL TENSILE FAILURE STRAIN (E1T) (FIBER DIRECTION)
CompressiveStrain1=-0.011			!LONGITUDINAL COMPRESSIVE FAILURE STRAIN (E1C)
TensileStrain2=0.0081				!TRANSVERSE TENSILE FAILURE STRAIN (E2T)
CompressiveStrain2=-0.032			!TRANSVERSE COMPRESSIVE FAILURE STRAIN (E2C)
TensileStrain3=0.007				!THROUGH-THICKNESS TENSILE FAILURE STRAIN (E3T)
CompressiveStrain3=-0.032			!THROUGH-THICKNESS COMPRESSIVE FAILURE STRAIN (E3C)

ShearStrain12=0.05					!IN-PLANE SHEAR FAILURE STRAIN (Gama12)
ShearStrain13=0.05					!TRANSVERSE SHEAR FAILURE STRAIN (Gama13)
ShearStrain23=0.021					!THROUGH-THICKNESS SHEAR FAILURE STRAIN (Gama23)

Density=1600						!MATERIAL DENSITY

Theta1=0
Theta2=90
Theta3=45
Theta4=-45

MPTEMP,,,,,,,,   
MPTEMP,1,0   
MPDATA,EX,1,,ElasticModulus1 
MPDATA,EY,1,,ElasticModulus2 
MPDATA,EZ,1,,ElasticModulus3 
MPDATA,PRXY,1,,PoissonsRatio12  
MPDATA,PRYZ,1,,PoissonsRatio23 
MPDATA,PRXZ,1,,PoissonsRatio13 
MPDATA,GXY,1,,ShearModulus12
MPDATA,GYZ,1,,ShearModulus23
MPDATA,GXZ,1,,ShearModulus13
MPDATA,DENS,1,,Density

TB,FCLI,1,1,9,1
TBTEMP,0
TBDATA,1,TensileStrength1
TBDATA,2,CompressiveStrength1
TBDATA,3,TensileStrength2
TBDATA,4,CompressiveStrength2
TBDATA,5,TensileStrength3
TBDATA,6,CompressiveStrength3
TBDATA,7,ShearStrength12 
TBDATA,8,ShearStrength23
TBDATA,9,ShearStrength13

TB,FCLI,1,1,9,2
TBTEMP,0
TBDATA,1,TensileStrain1
TBDATA,2,CompressiveStrain1
TBDATA,3,TensileStrain2
TBDATA,4,CompressiveStrain2
TBDATA,5,TensileStrain3
TBDATA,6,CompressiveStrain3
TBDATA,7,ShearStrain12 
TBDATA,8,ShearStrain23
TBDATA,9,ShearStrain13

FC,1,S,XTEN,TensileStrength1
FC,1,S,YTEN,TensileStrength2
FC,1,S,ZTEN,TensileStrength3
FC,1,S,XCMP,CompressiveStrength1
FC,1,S,YCMP,CompressiveStrength2
FC,1,S,ZCMP,CompressiveStrength3
FC,1,S,XY,ShearStrength12
FC,1,S,YZ,ShearStrength23
FC,1,S,XZ,ShearStrength13

FC,1,EPEL,XTEN,TensileStrain1
FC,1,EPEL,YTEN,TensileStrain2
FC,1,EPEL,ZTEN,TensileStrain3
FC,1,EPEL,XCMP,CompressiveStrain1
FC,1,EPEL,YCMP,CompressiveStrain2
FC,1,EPEL,ZCMP,CompressiveStrain3
FC,1,EPEL,XY,ShearStrain12
FC,1,EPEL,YZ,ShearStrain23
FC,1,EPEL,XZ,ShearStrain13


!*******************************************************************************************************
!								BEAM THICKNESS AND PLY NUMBER PROPERTIES
!*******************************************************************************************************

ThicknessTop=(LaminatePlyNumberTop*ThicknessPly)*12				!TOP THICKNESS (1 mm)*n
ThicknessWeb=(LaminatePlyNumberWeb*ThicknessPly)*4				!WEB THICKNESS (0,5 mm)*n
ThicknessBottom=(LaminatePlyNumberBottom*ThicknessPly)*12		!BOTTOM THICKNESS (1 mm)*n


ThicknessFlangeTotal=ThicknessTop+ThicknessBottom				!TOTAL FLANGE THICKNESS

PlyNumberTop=ThicknessTop/ThicknessPly							!TOTAL PLY NUMBER OF THE TOP FLANGE
PlyNumberWeb=ThicknessWeb/ThicknessPly							!TOTAL PLY NUMBER OF THE WEB
PlyNumberBottom=ThicknessBottom/ThicknessPly					!TOTAL PLY NUMBER OF THE BOTTOM FLANGE

!*******************************************************************************************************
!										SHELL THICKNESSES
!*******************************************************************************************************

SECT,1,shell,,Top Flange 
*DO,ply,1,PlyNumberTop/LaminatePlyNumberTop,1
SECDATA, ThicknessPly,1,Theta3,3    
SECDATA, ThicknessPly,1,Theta1,3   
SECDATA, ThicknessPly,1,Theta4,3 
SECDATA, ThicknessPly,1,Theta1,3    
SECDATA, ThicknessPly,1,Theta1,3    
SECDATA, ThicknessPly,1,Theta4,3 
SECDATA, ThicknessPly,1,Theta1,3   
SECDATA, ThicknessPly,1,Theta3,3    
SECOFFSET,TOP    
SECCONTROL,,,, , , , 
*ENDDO

SECT,2,shell,,Web
*DO,ply,1,PlyNumberWeb/LaminatePlyNumberWeb,1
SECDATA, ThicknessPly,1,Theta3,3    
SECDATA, ThicknessPly,1,Theta4,3   
SECDATA, ThicknessPly,1,Theta4,3   
SECDATA, ThicknessPly,1,Theta3,3    
SECOFFSET,MID    
SECCONTROL,,,, , , ,
*ENDDO

SECT,3,shell,,Bottom Flange  
*DO,ply,1,PlyNumberBottom/LaminatePlyNumberBottom,1
SECDATA, ThicknessPly,1,Theta3,3    
SECDATA, ThicknessPly,1,Theta1,3   
SECDATA, ThicknessPly,1,Theta4,3 
SECDATA, ThicknessPly,1,Theta1,3 
SECDATA, ThicknessPly,1,Theta1,3    
SECDATA, ThicknessPly,1,Theta4,3 
SECDATA, ThicknessPly,1,Theta1,3   
SECDATA, ThicknessPly,1,Theta3,3    
SECOFFSET,BOT    
SECCONTROL,,,, , , , 
*ENDDO

*END

!*******************************************************************************************************
!										STARTING PARAMETERS
!*******************************************************************************************************

*CREATE,StartingPar

n=4        					!DIMENSION OF THE PROBLEM (Number of Variables)
Lt=3*n        			    !MINIMUM LENGTH OF A MARKOW CHAIN
LtdOld=1.5*Lt  				!THIS VALUE IS USED ONLY IN THE FIRST DECREMENT OF THE CONTROL PARAMETER 
LtdOldP=LtdOld
NN=9*n     					!NUMBER OF CONFIGURATIONS

*DIM,Acost,ARRAY,NN   		!COST VALUE ARRAY, NN IS THE DIMENSION 
*DIM,order,ARRAY,NN   		!ORDER ARRAY, NN IS THE DIMENSION 
*DIM,Ah_rib,ARRAY,NN  	 	!INITIAL RIB HEIGHT ARRAY, NN IS THE DIMENSION 
*DIM,Aw_rib,ARRAY,NN  	 	!INITIAL RIB WIDTH ARRAY, NN IS THE DIMENSION 
*DIM,AL_rib,ARRAY,NN   		!INITIAL RIB LENGTH ARRAY, NN IS THE DIMENSION 
*DIM,Atheta_rib,ARRAY,NN   	!INTIAL ANGLE OF RIB ARRAY, NN IS THE DIMENSION
*DIM,ATsiWu,ARRAY,NN


ckini=200      				!INITIAL TEMPERATURE MODIFY FOR EVERY DIFFERENT MODEL
ck=ckini
alfamin=0.8  				!MINIMUM VALUE OF THE DECREMENT OF THE CONTROL PARAMETER
alfamax=0.9999  			!MAXIMUM VALUE OF THE DECREMENT OF THE CONTROL PARAMETER
epsin=0.00000002 			!STOPPING CRITERION (TEMPERATURE IS LOW ENOUGH)
epsinOut=0.0000002 			!STOPPING CRITERION (THE DIFFERENCE BETWEEN BEST AND WORST POINT IS NEGLIGIBLE)
NoIt=0         				!TOTAL NUMBER OF ITERATIONS EXCEPT GEOMETRICALLY UNACCEPTABLE MOVES
itotal=0       				!TOTAL NUMBER OF ITERATIONS 
mc=0           				!MARKOW CHAIN NUMBER
ann=-30        				!PARAMETER TO HANDLE VERY LOW NUMBERS

MassFirst=0					!INITIAL MASS
ObjectiveInitial=0			!INITIAL OBJECTIVE FUNCTION VALUE

!Defining Constraint Box for the Limits of the Rib
y4_min=-430e-3
y4_max=-20e-3
z4_min=-350e-3
z4_max=-22e-3

Hrib_min=0.005					!MINIMUM INITIAL NUMBER OF RIBS
Hrib_max=0.025					!MAXIMUM INITIAL NUMBER OF RIBS
Wrib_min=40e-3	 				!MINIMUM THICKNESS (WIDTH)
Wrib_max=100e-3     			!MAXIMUM THICKNESS (WIDTH)
L_rib_min=100e-3				!MINIMUM LENGTH
L_rib_max=300e-3				!MAXIMUM LENGTH 
Theta_rib_min=30				!MINIMUM ANGLE
Theta_rib_max=90				!MAXIMUM ANGLE

offset=-11e-4

rshrib_ini=(Hrib_max-Hrib_min)/3     	
rswrib_ini=(Wrib_max-Wrib_min)/3		
rsL_rib_ini=(L_rib_max-L_rib_min)/3		
rsthetarib_ini=(Theta_rib_max-Theta_rib_min)/3

rshrib=rshrib_ini			
rswrib=rswrib_ini			
rslrib=rsL_rib_ini			
rsthetarib=rsthetarib_ini

*END

!*******************************************************************************************************
!											INITIAL SHAPE
!*******************************************************************************************************

*CREATE,initialKey

!*******************************************************************************************************
!										SHAPE PROPERTIES OF THE BEAM
!*******************************************************************************************************

tc=2e-3			!thickness of web
ts=12e-3		!thickness of flange
H=500e-3		!Height of cross section 
Length=2000e-3	!Length of Beam
Wf=200e-3  		!Width of flange
LS=Length/10

!*******************************************************************************************************
!										SHAPE PROPERTIES OF THE RIB
!*******************************************************************************************************

H_rib=(Hrib_min+Hrib_max)/2			
W_rib=(Wrib_min+Wrib_max)/2			
L_rib=(L_rib_min+L_rib_max)/2			
Theta_rib=(Theta_rib_min+Theta_rib_max)/2

*END


!*******************************************************************************************************
!						DEFINITION OF THE KEYPOINTS AND GENERATION GEOMETRY
!*******************************************************************************************************

*CREATE,Definitionkey
/prep7

!*******************************************************************************************************
!											WEB AREA
!*******************************************************************************************************

L=Length
L1=Wf
L2=LS

K,1,-L1/2,0,(-L/2-L2)
K,2,0,0,(-L/2-L2)
K,3,L1/2,0,(-L/2-L2)
K,4,0,(-H+2*ts),(-L/2-L2)
K,5,-L1/2,(-H+2*ts),(-L/2-L2)
K,6,L1/2,(-H+2*ts),(-L/2-L2)

K,7,-L1/2,0,-L/2
K,8,0,0,-L/2
K,9,L1/2,0,-L/2
K,10,0,(-H+2*ts),-L/2
K,11,-L1/2,(-H+2*ts),-L/2
K,12,L1/2,(-H+2*ts),-L/2

K,13,-L1/2,0,0
K,14,0,0,0
K,15,L1/2,0,0
K,16,0,(-H+2*ts),0
K,17,-L1/2,(-H+2*ts),0
K,18,L1/2,(-H+2*ts),0

K,19,-L1/2,0,L/2
K,20,0,0,L/2
K,21,L1/2,0,L/2
K,22,0,(-H+2*ts),L/2
K,23,-L1/2,(-H+2*ts),L/2
K,24,L1/2,(-H+2*ts),L/2

K,25,-L1/2,0,(L/2+L2)
K,26,0,0,(L/2+L2)
K,27,L1/2,0,(L/2+L2)
K,28,0,(-H+2*ts),(L/2+L2)
K,29,-L1/2,(-H+2*ts),(L/2+L2)
K,30,L1/2,(-H+2*ts),(L/2+L2)

K,100,0,(-(H-2*ts)/2),0

*END

!*******************************************************************************************************
!					GENERATION OF THE AREAS AND CREATION OF THE LOCAL COORDINATE SYSTEMS
!*******************************************************************************************************

*CREATE,Areas

!---------------------------------------------------------------------------------------------------------
!									LOCAL COORDINATE SYSTEMS
!---------------------------------------------------------------------------------------------------------


*get,kp_max,KP,0,NUM,MAX
k,kp_max+1,0,0,-L/2
k,kp_max+2,0,0,(-L/2+1)
k,kp_max+3,L1/2,0,-L/2
k,kp_max+4,0,-(H-2*ts),-L/2
k,kp_max+5,0,-(H-2*ts),(-L/2+1)
k,kp_max+6,L1/2,-(H-2*ts),-L/2
k,kp_max+7,0,-(H-2*ts)/2,-L/2
k,kp_max+8,0,-(H-2*ts)/2,(-L/2+1)
k,kp_max+9,0,0,-L/2
k,kp_max+10,0,-(H-2*ts),L/2
k,kp_max+11,L1/2,-(H-2*ts),L/2
k,kp_max+12,0,-(H-2*ts)/2,L/2


CSKP,11,0,kp_max+1,kp_max+2,kp_max+3,    	!top flange coordinate
CSKP,12,0,kp_max+4,kp_max+5,kp_max+6, 		! bottom flange coordinate
CSKP,13,0,kp_max+7,kp_max+8,kp_max+9,    	! web coordinate
CSKP,14,0,kp_max+4,kp_max+6,kp_max+7,    	!edge at x=-l/2 
CSKP,15,0,kp_max+10,kp_max+11,kp_max+12,    !edge at x=L/2

!*******************************************************************************************************
!												WEB AREA
!*******************************************************************************************************

!Define Areas by Using Keyponts
csys,0

A,2,8,10,4
A,8,14,100,16,10
A,14,20,22,16,100
A,20,26,28,22

A,1,2,8,7
A,7,8,14,13
A,13,14,20,19
A,19,20,26,25
A,2,3,9,8
A,8,9,15,14
A,14,15,21,20
A,20,21,27,26

A,5,4,10,11
A,11,10,16,17
A,17,16,22,23
A,23,22,28,29
A,4,6,12,10
A,10,12,18,16
A,16,18,24,22
A,22,24,30,28

Asum
*get,Area0,AREA,0,AREA
!*******************************************************************************************************
!												Generate Ribs
!*******************************************************************************************************
csys,0
WPCSYS,-1,0 

rib1=0.20
rib2=0.38

k,151,offset,-(H-2*ts)/2+L_rib/2*sin(theta_rib)+W_rib/2*cos(theta_rib),L_rib/2*cos(theta_rib)-W_rib/2*sin(theta_rib)-rib1
k,152,H_rib,-(H-2*ts)/2+L_rib/2*sin(theta_rib),L_rib/2*cos(theta_rib)-rib1
k,153,offset,-(H-2*ts)/2+L_rib/2*sin(theta_rib)-W_rib/2*cos(theta_rib),L_rib/2*cos(theta_rib)+W_rib/2*sin(theta_rib)-rib1 

spline,151,152,153,,,,,1,-1,  ,-1,1

k,154,H_rib,-(H-2*ts)/2-L_rib/2*sin(theta_rib),-L_rib/2*cos(theta_rib)-rib1
l,152,154
adrag,51,52,,,,,53
    myerror1=_status
    *if,myerror1,ge,0.99,then
      reborder='yes'
    *endif
Ldele,53

k,159,offset,-(H-2*ts)/2+1.3*(L_rib/2+W_rib/2)*sin(theta_rib),1.3*(L_rib/2+W_rib/2)*cos(theta_rib)-rib1 
cskp,101,0,33,153,151,1,1
csys,101
kdele,154

spline,152,159,,,,,-1,0,0
adrag,51,52,,,,,53
    myerror2=_status
    *if,myerror2,ge,0.99,then
      reborder='yes'
    *endif

csys,0
k,162,offset,-(H-2*ts)/2-1.3*(L_rib/2+W_rib/2)*sin(theta_rib),-1.3*(L_rib/2+W_rib/2)*cos(theta_rib)-rib1
k,164,0,-(H-2*ts)/2-1.3*(L_rib/2+W_rib/2)*sin(theta_rib),-1.3*(L_rib/2+W_rib/2)*cos(theta_rib)-rib1
k,163,0,-(H-2*ts)/2+1.3*L_rib/2*sin(theta_rib),1.3*L_rib/2*cos(theta_rib)-rib1

csys,101
spline,32,162,,,,,1,0,0
adrag,54,57,,,,,64
    myerror3=_status
    *if,myerror3,ge,0.99,then
      reborder='yes'
    *endif

Ldele,64

KWPLAN,-1,     164,     163,     8
    myerror4=_status
    *if,myerror4,ge,0.99,then
      reborder='yes'
    *endif
	
block,-1.2*W_rib,1.6*L_rib+1.2*W_rib,-1.2*W_rib,1.2*W_rib,0,1.2*W_rib

asel,s,area,,21,26
asbv,all,1,,delete,delete
    myerror5=_status
    *if,myerror5,ge,0.99,then
      reborder='yes'
    *endif

allsel,all
aplot

al,88,89,92,85,83,91
    myerror6=_status
    *if,myerror6,ge,0.99,then
      reborder='yes'
    *endif
asba,2,21,,delete,delete
    myerror7=_status
    *if,myerror7,ge,0.99,then
      reborder='yes'
    *endif
allsel,all


csys,0
WPCSYS,-1,0 



agen,3,33,38,,,,-rib2,,0
    myerror8=_status
    *if,myerror8,ge,0.99,then
      reborder='yes'
    *endif
al,60,62,65,57,54,64
asba,22,39,,delete,delete
    myerror9=_status
    *if,myerror9,ge,0.99,then
      reborder='yes'
    *endif
al,73,75,78,70,68,77
asba,40,22,,delete,keep
asba,1,22,,delete,delete
    myerror10=_status
    *if,myerror10,ge,0.99,then
      reborder='yes'
    *endif
allsel,all

    myerror11=_status
    *if,myerror11,ge,0.99,then
      reborder='yes'
    *endif

FLST,3,18,5,ORDE,4  
FITEM,3,2   
FITEM,3,21  
FITEM,3,23  
FITEM,3,-38 
ARSYM,Z,P51X, , , ,0,0 

    myerror12=_status
    *if,myerror12,ge,0.99,then
      reborder='yes'
    *endif


al,131,133,136,128,126,135
asba,3,57,,delete,delete
    myerror13=_status
    *if,myerror13,ge,0.99,then
      reborder='yes'
    *endif
al,105,107,110,102,100,109
asba,58,3,,delete,delete
    myerror14=_status
    *if,myerror14,ge,0.99,then
      reborder='yes'
    *endif
al,118,120,123,115,113,122
asba,57,3,,delete,keep
asba,4,3,,delete,delete
    myerror15=_status
    *if,myerror15,ge,0.99,then
      reborder='yes'
    *endif
allsel,all

aglue,all

    myerror17=_status
    *if,myerror17,ge,0.99,then
      reborder='yes'
    *endif
/Output,temp1,out,,Append
*VWRITE,'17  ',myerror17
(A4,F5.0)              
*VWRITE,jk,H_rib,W_rib,L_rib,Theta_rib,
(F5.0,F14.9,F14.9,F14.9,F14.9)  
 /Output         
*Cfclos

*END


!==============================================================!=========================================
!			GENERATION OF THE MESHES, APPLYING THE LOAD AND STABILIZATION BOUNDARY CONDITION
!==============================================================!=========================================
*CREATE,Meshes
/prep7
!---------------------------------------------------------------------------------------------------------
!		MESH
!---------------------------------------------------------------------------------------------------------

*do,i,5,12,1
CM,_Y,AREA  
ASEL, , , ,       i
CM,_Y1,AREA 
CMSEL,S,_Y  
!*  
CMSEL,S,_Y1 
AATT,       1, ,   1,       11,   1  
CMSEL,S,_Y  
CMDELE,_Y   
CMDELE,_Y1  
*enddo

!--------------------------------------------------------------------------------------------------------- 

*do,i,13,20,1
CM,_Y,AREA  
ASEL, , , ,       i
CM,_Y1,AREA 
CMSEL,S,_Y  
!*  
CMSEL,S,_Y1 
AATT,       1, ,   1,       12,   3  
CMSEL,S,_Y  
CMDELE,_Y   
CMDELE,_Y1  
*enddo

!--------------------------------------------------------------------------------------------------------- 

*do,i,1,4,1
CM,_Y,AREA  
ASEL, , , ,       i
CM,_Y1,AREA 
CMSEL,S,_Y  
!*  
CMSEL,S,_Y1 
AATT,       1, ,   1,       13,   2  
CMSEL,S,_Y  
CMDELE,_Y   
CMDELE,_Y1  
*enddo  

!---------------------------------------------------------------------------------------------------------
numcmp,area
*get,AreaTotalNo,area,0,count
*do,i,21,AreaTotalNo,1
CM,_Y,AREA  
ASEL, , , ,       i
CM,_Y1,AREA 
CMSEL,S,_Y  
!*  
CMSEL,S,_Y1 
AATT,       1, ,   1,       13,   2  
CMSEL,S,_Y  
CMDELE,_Y   
CMDELE,_Y1  
*enddo


!---------------------------------------------------------------------------------------------------------
!MESHING
!---------------------------------------------------------------------------------------------------------

Mesh_size=8e-3
Esize,8e-3				!Specifies the Default Number of Line Divisions
Amesh,All				!Generate Nodes and Area Elements within All Area by Using Defines Size in Esize Command
    myerror18=_status
    *if,myerror18,ge,2.99,then
      reborder='yes'
    *endif
/Output,temp1,out,,Append
*VWRITE,'18  ',myerror18
(A4,F5.0)              
/Output
*Cfclos

Asum
*get,TotalArea,AREA,0,AREA

*END


!*******************************************************************************************************
!											SOLUTION
!*******************************************************************************************************
*CREATE,Solution
/prep7
!*******************************************************************************************************
!								APPLIED LOAD AND BOUNDARY CONDITIONS
!*******************************************************************************************************

csys,0
!Defines DOF Constraints on Lines

DL,34,,UX,0
DL,34,,UY,0
DL,44,,UX,0
DL,44,,UY,0
DL,38,,UX,0
DL,38,,UY,0
DL,38,,UZ,0
DL,48,,UX,0
DL,48,,UY,0
DL,48,,UZ,0

DK,7,UX,0
DK,9,UX,0
DK,21,UX,0
DK,19,UX,0
DK,7,ROTZ,0
DK,9,ROTZ,0
DK,21,ROTZ,0
DK,19,ROTZ,0
!

!---------------------------------------------------------------------------------------------------------
!       				Master node for CERIG is created
!---------------------------------------------------------------------------------------------------------
K,5e6,0,0,0
ET,3,MASS21
R,1,0,0,0,0,0,0,
NUMSTR,NODE,5e6
KMESH,       5e6

!---------------------------------------------------------------------------------------------------------
!       				Rigid Region is created
!---------------------------------------------------------------------------------------------------------
NSEL,S,LOC,Z,(-(Mesh_size)/2),((Mesh_size)/2)
NSEL,R,LOC,Y,0,ts
NSEL,U, , ,    5e6 
CM,load_nodes,NODE,
CMSEL,S,load_nodes
NSEL,A, , ,    5e6 
CERIG,5e6,ALL,UY,RXYZ

ALLSEL

Force=-1   

F,5e6,FY,Force
Alls

 

FINISH  
/Solu
Alls
Antype,0							!Static Analysis
Solve
    myerror19=_status
    *if,myerror19,ge,1.99,then
      reborder='yes'
    *endif
/Output,temp1,out,,Append
*VWRITE,'19  ',myerror19
(A4,F5.0)              
/Output
*Cfclos

Finish

/solu
ANTYPE,1							!Linear Buckling Analysis
!*  
BUCOPT,LANB,6,0,0,CENTER
Pstres,on   
/STATUS,SOLU
SOLVE   
    myerror20=_status
    *if,myerror20,ge,1.99,then
      reborder='yes'
    *endif
/Output,temp1,out,,Append
*VWRITE,'20  ',myerror20
(A4,F5.0)              
/Output
*Cfclos

FINISH  

SAVE,'buckling','db'
FINISH

*END


!*******************************************************************************************************
!								CREATING THE MASS AND OBJECTIVE FUNCTIONS
!*******************************************************************************************************

*CREATE,costweightobj,mac


FINISH

/post1

!--------------------------------------------------------------------------------------------------------- 
!Mass0  =  Mass of I-beam without RIB
!Area0  =  Area of I-beam without RIB
!TotalArea = Total Area of I-beam With RIB

MassFirst=MassFirst+1
!Flange_AREA=L1*(Length+LS)*2
!Web_AREA=(H-2*ts)*(Length+LS)
!Rib_AREA=H_rib*((0.25*3.14*W_rib*W_rib)+(W_rib*L2))

Mass0=Density*9.81*Area0

!--------------------------------------------------------------------------------------------------------- 
!BuckleLF0 =  Buckling Load Factor of I-beam Without Rib
!BuckleLF  =  Buckling Load Factor of I-beam With Rib

BuckleLF0=5721.91

*GET,BuckleLF,MODE,1,FREQ

Mass=TotalArea*9.81*Density

*if,BuckleLF,le,0.0,then
   reborder='yes'
/Output,temp1,out,,Append
*VWRITE                                         
(T2,'Negative buckling factor',)                    
/Output
*Cfclos
*endif

ObjectiveInitial=ObjectiveInitial+1
Obj_F=0.85*(BuckleLF0/BuckleLF) + 0.15*(Mass/Mass0)

!---------------------------------------------------------------------------------------------------------

!FCTYP,ADD,TWSI
!nsort,fail,TWSI,0,0
!*get,MaxFail,sort,0,max

!---------------------------------------------------------------------------------------------------------
/OUTPUT,Objective_out,out,,APPEND
*VWRITE,jk,Obj_F,Mass0,Mass,BuckleLF,H_rib,W_rib,L_rib,Theta_rib,MaxFailTsaiWu
(F5.0,F14.6,F14.6,F14.6,F15.4,F12.7,F12.7,F12.7,F12.7,F9.5)
/OUTPUT
*CFCLOS

	/GRAPHICS,POWER
FINISH


*END

!*******************************************************************************************************
!					CALCULATING SHAPE PROPERTIES AND COST VALUES FOR THE CURRENT SHAPE
!*******************************************************************************************************

*CREATE,RefreshC

H_rib=Ah_rib(CurrentShape)			!INITIAL WAVE LENGTH
W_rib=Aw_rib(CurrentShape)			!INITIAL AMPLITUDE
L_rib=AL_rib(CurrentShape)			!FINAL AMPLITUDE
Theta_rib=Atheta_rib(CurrentShape)

*END

!*******************************************************************************************************
!							GENERATION OF RANDOM SHAPE PROPERTIES
!*******************************************************************************************************

*CREATE,RanPoint,mac

*DO,wwOut,1,10000

	*DO,ww,1,10000
		Step=rshrib*RAND(-1,1)
		TempH_rib=H_rib+Step
		*IF,TempH_rib,GT,Hrib_max,CYCLE
		*IF,TempH_rib,LT,Hrib_min,CYCLE
		*EXIT
	*ENDDO
	H_rib=TempH_rib

	!*IF,ww,GT,9998,THEN
	!	reborder='yes'
	!*ENDIF

	*DO,ww,1,10000
		Step=rswrib*RAND(-1,1)
		TempW_rib=W_rib+Step
		*IF,TempW_rib,GT,Wrib_max,CYCLE
		*IF,TempW_rib,LT,Wrib_min,CYCLE
		*EXIT
	*ENDDO
	W_rib=TempW_rib

	!*IF,ww,GT,9998,THEN
	!	reborder='yes'
	!*ENDIF
	
	*DO,ww,1,10000
		Step=rslrib*RAND(-1,1)
		TempL_rib=L_rib+Step
		*IF,TempL_rib,GT,L_rib_max,CYCLE
		*IF,TempL_rib,LT,L_rib_min,CYCLE
		*EXIT
	*ENDDO
	L_rib=TempL_rib

	!*IF,ww,GT,9998,THEN
	!	reborder='yes'
	!*ENDIF

	*DO,ww,1,10000
		Step=rsthetarib*RAND(-1,1)
		TempTheta_rib=Theta_rib+Step
		*IF,TempTheta_rib,GT,Theta_rib_max,CYCLE
		*IF,TempTheta_rib,LT,Theta_rib_min,CYCLE
		*EXIT
	*ENDDO
	Theta_rib=TempTheta_rib

	!*IF,ww,GT,9998,THEN
	!	reborder='yes'
	!*ENDIF
	
	
	Y4_rib=-(H-2*ts)/2+1.3*(L_rib/2+W_rib/2)*sin(Theta_rib)
	Z4_rib=1.3*(L_rib/2+W_rib/2)*cos(Theta_rib)-0.2


		*IF,Y4_rib,GT,y4_max,CYCLE
		*IF,Y4_rib,LT,y4_min,CYCLE
		*IF,Z4_rib,GT,z4_max,CYCLE
		*IF,Z4_rib,LT,z4_min,CYCLE
		*EXIT
	*ENDDO
	

	*IF,ww,GT,9998,THEN
		reborder='yes'
	*ENDIF

*END


!*******************************************************************************************************
!								FINDING THE BEST AND WORST POINTS
!*******************************************************************************************************

*CREATE,Reorder

order(1)=1
*DO,jj,2,NN
	order(jj)=jj
	*DO,ii,1,jj-1
		*IF,Acost(jj),LT,Acost(order(ii)),THEN
			*DO,kk,jj,ii+1,-1
				order(kk)=order(kk-1)
			*ENDDO
			order(ii)=jj
			*EXIT
		*ENDIF
	*ENDDO
*ENDDO
Best=order(1)
Worst=order(NN)
Worse=order(NN-3*n)
WorstBest=order(NN-3*n+1)

*END


!---------------------------------------------------------------------------------------------------------

!---------------------------------------------------------------------------------------------------------

!---------------------------------------------------------------------------------------------------------


!START
!/clear,start 					!if you are restarting erase the above lines
*DO,rknpp,1,1803
	Temp=RAND(0,1) 					!These lines prevent generating the same random moves 
*ENDDO 								!(ANSYS seems to Generate the same random points every time the program starts)
/input,,parm          			!if you are restarting activate this line
repeatNo=mod(mc,20)+1			!this parameter is calculated to obtain a better randomization
s=1
parsav,all                	    !thus at each restart, different random numbers are generated


!*******************************************************************************************************
!											Outer loop (Markov Chain)
!*******************************************************************************************************

*DOWHILE,s
FINISH
/CLEAR,START
/INPUT,,parm

*if,mod(mc,50),eq,0,then
/Output,Param,out,,Append
*VWRITE                                         
(T2,'mc   Difcost MarkL MrkLA AccMv BestBuckF BestObj BestMass  Worse    worst   StepSize Temperature best Noim RatioAcc RatioSt')            
/Output
*Cfclos
*endif

	*do,rkn,1,repeatNo*50   
		temp=rand(0,1)   	!these lines prevent generating the same random moves
	*enddo          	!(ANSYS seems to generate the same random points every time the program starts) 

*if,exitLoop,eq,'true',then  

mc=mc+1
exitOut='false'
exitLoop='false'
Noim=0          !No of imrovements made on the worse point
AccMov=0        !No of accepted moves   
in=0            !counter of the inner loop (Markow chain) iterations

LtdNew=nint(Lt*(1+2*(rshrib/rshrib_ini)))     !current Markow length
LtdNewP=LtdNew     							  !if no point is found in the Morkow chain which is better than the best point
*ENDIF

!*******************************************************************************************************
!											Inner loop (Markov Chain)
!*******************************************************************************************************

ss=1
*DOWHILE,ss
in=in+1
NoIt=NoIt+1

!*****GENERATING A NEW POINT*****
sss=1
parsav,all
*DOWHILE,sss
  FINISH

  /clear,start  
  /input,,parm  
  *use,Initial   				!initialization of the program 

  /prep7
  itotal=itotal+1
  reborder='no'
  *do,rkn,1,repeatNo  		 !these lines prevent generating the same random moves
	temp=rand(0,1)     		 !(ANSYS seems to generate the same random points every time the program starts)
  *enddo 

  CurrentShape=nint(rand(0.51,NN+0.49))  !choosing randomly the current point in A
  *use,refreshC                          !calculating x,y and cost values for the current point
  RanPoint  							   !Generating a random point

  *USE,Definitionkey
  *IF,reborder,EQ,'no',THEN
    *use,Areas  
  *endif
  *IF,reborder,EQ,'no',THEN
    *use,Meshes 
  *endif
  *IF,reborder,EQ,'no',THEN
	*use,solution
  *endif
  *IF,reborder,EQ,'no',THEN
    costweightobj   					   !calculating the cost of the new point 
   *endif
  *IF,reborder,EQ,'no',THEN
    sss=0
  *endif
*enddo

exponent=(Acost(worstBest)-Obj_F)/ck
*if,Obj_F,le,Acost(worstBest),then
	AcceProb=1
*elseif,exponent,le,ann,then
	AcceProb=exp(ann)
*else
	AcceProb=exp(exponent)
*endif

randomNo=rand(0,1)
*if,randomNo,le,AcceProb,then   		 !Accept the move
  AccMov=AccMov+1    
  *if,Obj_F,lt,Acost(best),then
    LtdNewP=in
    exitLoop='true'  					 !Since a point better than the best one if found, the the current loop is ended
	BestObj_F=Obj_F
	BestTsaiWuIndexMax = MaxFail
	BestBuckleLF=BuckleLF
	BestMass=Mass	
	Best_H = H_rib
	Best_w = W_rib  
	Best_L = L_rib
	Best_theta=theta_rib
/Output,BestConfigurations,out,,Append
*VWRITE,mc,AcostB,BestBuckleLF,BestMass,BestTsaiWuIndexMax,Best_H,Best_w,Best_L,Best_theta
(F5.0,F9.5,F11.2,F11.3,F9.6,F9.6,F9.6,F9.6,F10.6)
/Output
*Cfclos

  *endif
  *if,Obj_F,lt,Acost(worse),then  	!This means that the configuration better than the worst but worse than the rest has been improved
    Noim=Noim+1                    
  *endif
  mm=NN-nint(rand(0.501,n+0.499))+1
  worst=order(mm)
	Acost(worst)=Obj_F
    Ah_rib(worst)=H_rib
    Aw_rib(worst)=W_rib  
    AL_rib(worst)=L_rib
	Atheta_rib(worst)=theta_rib
	ATsiWu(worst)=MaxFailTsaiWu
  *use,Reorder   !Finding the best and worst points
*endif

*if,in,ge,LtdNew-1E-6,then
  exitLoop='true'
*endif

*if,exitLoop,eq,'true',then
  ss=0
*endif
*enddo

*if,Noim,lt,0.1*LtdNewp,then
  rshrib=0.9*rshrib   
  rswrib=0.9*rswrib
  rslrib=0.9*rslrib
  rsthetarib=0.9*rsthetarib
*endif         


!*******************************************************************************************************
!											End of inner loop (Markov chain)
!*******************************************************************************************************

!Decrement of the temperature

RatioAcc=AccMov/in
RatioStep=(rshrib/rshrib_ini)+0.01
*if,RatioAcc,gt,RatioStep,then
  alfa=alfamin
*else
 alfa=alfamax
*endif
*if,LtdNewP,ne,LtdNew,then
  alfa=1
*endif

Difcost=Acost(worse)-Acost(best)
AcostWB=Acost(worse)
AcostB=Acost(best)
AcostWst=Acost(worst)

/Output,Param,out,,Append
*VWRITE,mc,Difcost,LtdNew,LtdNewP,AccMov,BestBuckleLF,AcostB,BestMass,AcostWB,AcostWst,rshrib,ck,best,Noim,RatioAcc,RatioStep
(F5.0,F9.6,F5.0,F5.0,F6.0,F10.2,F9.5,F10.2,F9.6,F9.6,F9.6,F12.4,F5.0,F5.0,F8.5,F8.5)
/Output
*Cfclos

!Control SA stopping criterion

*if,ck,lt,epsin,then
  *if,(AcostWst-AcostB),lt,epsinout,then
    exitOut='true'
  *endif
*endif

!Decrement of the step size

ck=alfa*ck

parsav,all
/copy,,parm,,mark%mc%,,,

*if,exitOut,eq,'true',then
s=0
*endif
*enddo

!Regenerating the best solution

/prep7
		ALLSEL,ALL
		ACLEAR,ALL
		ADELE,ALL
		LDELE,ALL
		KDELE,ALL
		EDELE,ALL      
    
/replot
/prep7
*use,Initial
*use,Definitionkey
CurrentShape=best
*use,refreshC    
*use,Areas  
*use,Meshes 
*use,solution 
costweightobj              
/post1
/graphics,full