import string
import os
#
#########################################################################
#
#  THIS SCRIPT RUNS A SEQUENCE OF PARAMETERIZED INPUT FILES TO STUDY THE
#  POSTBUCKLING BEHAVIOR OF A LINEAR ELASTIC, AXIALLY LOADED, 
#  CYLINDRICAL SHELL.  
#
#  PARAMETERS USED IN STUDY:
#
#  SHELL THICKNESS:                        thickness    (set by xthickness)
#  LENGTH OF THE CYLINDER:                 length       (set by xlength)
#  MEAN RADIUS OF THE CYLINDER:            radius       (set by xradius)
#  NUMBER OF NODES AROUND CIRCUMF:         node_circum  (set by nel_c)
#  NUMBER OF NODES ALONG LENGTH:           node_length  (set by nel_l)
#  POISSON'S RATIO:                        poisson      (set by xpoisson)
#  YOUNG'S MODULUS:                        young        (set by xyoung)
#  APPLIED LOAD (BUCKLING ANALYSIS):       tot_load     (set by xload)
#  INTERNAL PRESSURE:                      int_press    (set by press)
#  NUMBER OF BUCKLING MODES:               num_modes    (set by nmodes)
#
#  ADDITIONAL PARAMETERS FOR POST-BUCKLING ANALYSIS
#
#  APPLIED AXIAL LOAD (BASED ON            
#  LOWEST BUCKLING MODE):                  eig1_load    (set by xload1)
#  BUCKLING ANALYSIS RESULTS FILE NAME:    buckle_file  (set by bklfnames) 
#  *INCLUDE FILE WITH *IMPERFECTION DATA : imperf_file  (set by impfnames)
#
#  ADDITIONAL VARIABLES
#
#  IMPERFECTION SCALE FACTOR:              RadialImperfFactor
#  PRECRIBED RADIAL IMPERFECTION:          rad_imp
#  EIGENMODES TO SEED THE IMPERFECTION:    eigmodes
#  SCALE FACTOR ASSOC. WITH EACH MODE:     modefctr
#     -THIS IS A GUESS; THE TRUE FACTORS
#      WILL BE COMPUTED LATER.
#
#########################################################################
#
#  1. DEFINE THE SHELL PARAMETERS AND MAX-OUT-OF ROUND IMPERFECTION FACTOR

xthickness = 0.01
xlength = 2.0
xradius = 5.0
nel_c = 240
nel_l = 21
xpoisson = 0.3
xyoung = 30.0e6
xload = 1.0
press = 1.0
nmodes = 19

radialImperfFactor = 0.5

#  2. SPECIFY THE EIGENMODES WHICH WILL BE USED TO SEED THE IMPERFECTION
#     (TAKE ALL MODES, 1 THROUGH NMODES)

eigmodes = [ ]
for i in range(nmodes): eigmodes.append(i+1)

#  3. MISC. DEFINITIONS

zero = 0

#########################################################################
#
#       PERFORM EIGENVALUE BUCKLING ANALYSIS
#
#       CREATE THE STUDY

buckle = ParStudy(par=('thickness','length','radius',
                     'node_circum','node_length','poisson','young',
                     'tot_load','int_press','num_modes'))

#  DEFINE THE NAMES OF THE INPUT DECKS

names = ['cylsh_buck']

#  DEFINE THE PARAMETERS

buckle.define(DISCRETE,par='thickness')
buckle.define(DISCRETE,par='length')
buckle.define(DISCRETE,par='radius')
buckle.define(DISCRETE,par='node_circum')
buckle.define(DISCRETE,par='node_length')
buckle.define(DISCRETE,par='poisson')
buckle.define(DISCRETE,par='young')
buckle.define(DISCRETE,par='tot_load')
buckle.define(DISCRETE,par='num_modes')
buckle.define(DISCRETE,par='int_press')

#  SAMPLE THE PARAMETERS - INPUT THE APPROPRIATE VALUES

buckle.sample(VALUES,par='thickness',values=(xthickness))
buckle.sample(VALUES,par='length',values=(xlength))
buckle.sample(VALUES,par='radius',values=(xradius))
buckle.sample(VALUES,par='node_circum',values=(nel_c))
buckle.sample(VALUES,par='node_length',values=(nel_l))
buckle.sample(VALUES,par='poisson',values=(xpoisson))
buckle.sample(VALUES,par='young',values=(xyoung))
buckle.sample(VALUES,par='tot_load',values=(xload))
buckle.sample(VALUES,par='num_modes',values=(nmodes))
buckle.sample(VALUES,par='int_press',values=(press))

#  COMBINE THE SAMPLES INTO ANALYSES

buckle.combine(MESH,name='short')

for temp in names:

    #  GENERATE INPUT DECKS AND EXECUTION SCRIPT
    #  FOR VARIOUS TEMPLATES

    buckle.generate(template=temp)

    #  EXECUTE RUNS SEQUENTIALLY
    
    buckle.execute()

    #  GATHER RESULTS FOR FIRST BUCKLING MODE AND WRITE TO OUTPUT FILES
    
    buckle.gather(results='mode',variable='MODAL',step=2)
    buckle.report(PRINT,par=('length','radius','thickness',
                             'poisson','young'),results=('mode.2'))
                                 

    #   3. STORE LOWEST EIGENVALUE IN EIG_VALS; COMPUTE RAD_IMP.
    #   
    #   STORE THE RESULTS (.FIL) FILE NAMES OF EACH ANALYSIS (IN BKLFNAMES), 
    #   THE PARAMETER SETTINGS FOR EACH ANALYSIS (IN PARAMS), AND A 
    #   DESCRIPTIVE STRING OF CHARACTERS FOR EACH ANALYSIS (IN XNAMES). 
    #   STORE THE NUMBER OF NODES ALONG THE CIRCUMFERENCE (IN N_CIRC) AND 
    #   LENGTH (N_LGTH)

    eig_vals=[ ]        
    rad_imp =[ ]
    n_circ = [ ]
    n_lgth = [ ]
    bklfnames=[ ]
    xnames=[ ]
    params = []
    i = -1
    res = buckle.table.results
    for jname in buckle.job.keys():
        i = i+1
        des = buckle.job[jname].design

        thickness  = des[0]
        r_imp = radialImperfFactor*thickness
        length = des[1]
        radius = des[2]
        ncrc = des[3]
        nlen = des[4]
        value = res[i][0]
        xname = buckle.job[jname].designName
        root = buckle.job[jname].root
        bklfname = root + '_' + xname
        bklfnames.append(bklfname)
        xnames.append(xname)
        param = des
        params.append(param)
        
        eig_vals.append(value)
        rad_imp.append(r_imp)
        n_circ.append(ncrc)
        n_lgth.append(nlen)

#  INTERMEDIATE STEP:  RUN FORTRAN PROGRAM TO 
#  DETERMINE TRUE SCALE FOR IMPERFECTIONS
#
#  NOW CREATE THE INPUT FILE REQUIRED TO RUN THE FORTRAN PROGRAM;
#  INCLUDES THE FOLLOWING:
#
#  1. NAME ASSIGNED TO THE OUTPUT FILE CREATED BY THE PROGRAM
#     THIS FILE CONTAINS THE SCALE FACTOR FOR EACH EIGENMODE 
#     USED TO SEED THE IMPERFECTION AND WILL BE INCLUDED IN THE
#     POST BUCKLING ANALYSIS FILE.)
#  2. RESULTS FILE NAME FOR EIGENVALUE BUCKLING ANALYSIS
#  3. NUMBER OF CIRCUMF. AND LONG. NODES
#  4. PRESCRIBED RADIAL IMPERFECTION (FROM RAD_IMP)
#  5. LIST OF EIGENMODES (FROM EIGMODES) FOLLOWED BY ZERO
#  6. LIST OF "GUESS" SCALE FACTORS (FROM MODEFCTR)  
#

        file = 'max_round_input.dat'
        impfnames = [ ]
        names = [ ]
        i = -1
        for bklfname in bklfnames:
            i = i+1
            impfname = bklfname + str('_imp')
            impfnames.append(impfname)
            names.append(bklfname)
            
            modefctr = [ ]
            for m in eigmodes: modefctr.append(rad_imp[i])
            f1 = open(file,'w')
            f1.write(impfnames[i] + '\n')
            f1.write(names[i] + '\n')
            f1.write(str(n_circ[i]) + '\n')
            f1.write(str(n_lgth[i]) + '\n')
            f1.write(str(rad_imp[i]) + '\n')
            
            for imode in eigmodes: f1.write(str(imode) + '\n')
            f1.write(str(zero) + '\n')
            for ifctr in modefctr: f1.write(str(ifctr) + '\n')
            f1.close()

#  RUN THE FORTRAN PROGRAM
	    if os.name == 'nt':
                    os.system('cylsh_maximp.exe')
	    else:
                    os.system('./cylsh_maximp')

#  PERFORM THE POST-BUCKLING ANALYSIS
#
#  CREATE THE STUDY

riks = ParStudy(par=('thickness','length','radius',
                     'node_circum','node_length','poisson','young',
                     'int_press','eig1_load','buckle_file','imperf_file'))

#  DETERMINE THE APPROPRIATE LOAD LEVEL FOR THE ANALYSIS

xload1 = [ ]
for l in eig_vals:
    xl = float(xload)*float(l)
    xload1.append(xl)

#  DEFINE THE NAMES OF THE INPUT DECKS

names2 = ['cylsh_postbuck']

#  DEFINE THE PARAMETERS

riks.define(DISCRETE,par='thickness')
riks.define(DISCRETE,par='length')
riks.define(DISCRETE,par='radius')
riks.define(DISCRETE,par='node_circum')
riks.define(DISCRETE,par='node_length')
riks.define(DISCRETE,par='poisson')
riks.define(DISCRETE,par='young')
riks.define(DISCRETE,par='eig1_load')
riks.define(DISCRETE,par='buckle_file')
riks.define(DISCRETE,par='imperf_file')
riks.define(DISCRETE,par='int_press')

#  SAMPLE THE PARAMETERS - INPUT THE APPROPRIATE VALUES

for i in range(len(xload1)):

    riks.sample(VALUES,par='thickness',values=(params[i][0]))
    riks.sample(VALUES,par='length',values=(params[i][1]))
    riks.sample(VALUES,par='radius',values=(params[i][2]))
    riks.sample(VALUES,par='node_circum',values=(params[i][3]))
    riks.sample(VALUES,par='node_length',values=(params[i][4]))
    riks.sample(VALUES,par='poisson',values=(params[i][5]))
    riks.sample(VALUES,par='young',values=(params[i][6]))
    riks.sample(VALUES,par='int_press',values=(params[i][8]))
    riks.sample(VALUES,par='eig1_load',values=(xload1[i]))
    riks.sample(VALUES,par='buckle_file',values=(bklfnames[i]))
    riks.sample(VALUES,par='imperf_file',values=(impfnames[i]))

    #  COMBINE THE SAMPLES INTO ANALYSES
    #
    #  (APPEND XNAME TO THE JOBNAMES OTHERWISE THEY WILL NOT BE UNIQUE; 
    #   ALSO PROVIDES A MEANS OF ASSOCIATING THE POSTBUCKLING ANALYSIS 
    #   WITH ITS CORRESPONDING EIGENVALUE BUCKLING ANALYSIS.)

    riks.combine(MESH,name=xnames[i])

for temp in names2:

    #  GENERATE INPUT DECKS AND EXECUTION SCRIPT
    #  FOR VARIOUS TEMPLATES

    riks.generate(template=temp)

    #  EXECUTE RUNS SEQUENTIALLY

    riks.execute()
    riks.gather(results='lpf',variable='LPF',step=2,inc=LAST)

    #  GATHER RESULTS FOR ABAQUS QA PURPOSES ONLY

    riks.report(PRINT,par=('length','radius',
                           'thickness','poisson','young'),
                results=('lpf'))
    buck_res = buckle.table.results
    riks_res = riks.table.results
    resfile = open('cylsh_script.psr','w')
    for x in buck_res: resfile.write('%14.6g \n' % x[0])
    for y in riks_res: resfile.write('%14.2g \n' % y[0])