import string

#########################################################################
#
#  This script runs a sequence of parameterized input files to compare 
#  the relative accuracy of shell, elbow (with and without ovalization), 
#  and pipe elements. 
#
#  Parameters used in study:
#
#  wall thickness:                         wall_thick
#  prescribed displacements:               disp (set through dsp)
#  length of the straight section:         a
#  radius of the curved section:           bend_radius
#  outer radius of the structure:          outer_pipe_radius
#  number of elements in straight section: num_elem_s
#  number of elements in curved section:   num_elem_c
#  flexibility factor k:                   k_factor
#  Poissons ratio                          poisson
#  Young's modulus                         young
#  Output node                             n1 (set through nout)
#
#########################################################################

# SET PRESCRIBED DISPLACEMENT AND OUTPUT NODE
dsp  = 1.0
nout = 1

# DEFINE A FUNCTION TO OUTPUT DATA TO XYPLOT FILES
def customTable(file1, file2):
    for line in file1.readlines():
	print line
	nl = string.split(line,',')

	disp = float(nl[0])
	bend_radius = float(nl[1])
	wall_thick  = float(nl[2])
	outer_pipe_radius = float(nl[3])
	poisson = float(nl[4])
	rf = float(nl[6])

	mean_rad = outer_pipe_radius - wall_thick/2.0
	k = bend_radius*wall_thick/mean_rad**2/sqrt(1.e0 - poisson**2)
	k = 1.66e0/k

	outputstring = str(k) + ', ' + str(rf) + '\n'
	file2.write(outputstring)

# CREATE THE STUDY
pars = ('wall_thick','disp','a','bend_radius','outer_pipe_radius',
	'num_elem_s','num_elem_c','poisson','young','n1')
tp = ParStudy(par=pars, directory=OFF, verbose=ON, name='elbowtest')

# DEFINE THE NAMES OF THE INPUT DECKS
inputFiles=('elbowtest_shell','elbowtest_elbow6','elbowtest_pipek',
	    'elbowtest_elbow0','elbowtest_elbow3','elbowtest_pipe')

# DEFINE THE PARAMETERS
tp.define(DISCRETE, par=pars)

# SAMPLE BY VALUE
tp.sample(VALUES, par='wall_thick', values=(.08,.07,.06,.05,.04,.03,.02,.0125))
tp.sample(VALUES, par='disp', values=(dsp))
tp.sample(VALUES, par='a', values=(10.))
tp.sample(VALUES, par='bend_radius', values=(4.))
tp.sample(VALUES, par='outer_pipe_radius', values=(.5))
tp.sample(VALUES, par='num_elem_s', values=(25))
tp.sample(VALUES, par='num_elem_c', values=(25))
tp.sample(VALUES, par='poisson', values=(0.0))
tp.sample(VALUES, par='young', values=(28.1E6))
tp.sample(VALUES, par='n1', values=(nout))

# COMBINE THE SAMPLES INTO DESIGNS
tp.combine(MESH, name='open')

# RE-SAMPLE AND COMBINE
tp.sample(VALUES, par='disp', values=(-dsp))
tp.combine(MESH, name='close')

# BEGIN LOOP OVER INPUT FILES
for fileName in inputFiles:

    # GENERATE ANALYSIS JOB DATA
    # FOR VARIOUS TEMPLATES

    tp.generate(template=fileName)

    # EXECUTE ALL JOBS SEQUENTIALLY
    tp.execute(ALL)

    # PARAMETRIC STUDY OUTPUT AT END OF STEP 1
    tp.output(step=1, inc=LAST)

    # GATHER RESULTS FOR OUTPUT NODE AND WRITE TO OUTPUT FILES
    tp.gather(results='react', variable='RF', node=nout)

    fname_o = fileName + '_open_output'

    tp.report(PRINT, designSet='open',
	      par=('disp','bend_radius','wall_thick',
		   'outer_pipe_radius','poisson','young'),
	      file=fname_o, results=('react.1'))

    fname_c = fileName + '_close_output'

    tp.report(PRINT, designSet='close',
	      par=('disp','bend_radius','wall_thick',
		   'outer_pipe_radius','poisson','young'),
	      file=fname_c, results=('react.1'))

    tp.report(XYPLOT, designSet='open',
	      par=('disp','bend_radius','wall_thick',
		   'outer_pipe_radius','poisson','young'),
	      file=fname_o, results=('react.1'))

    tp.report(XYPLOT, designSet='close',
	      par=('disp','bend_radius','wall_thick',
		   'outer_pipe_radius','poisson','young'),
	      file=fname_c, results=('react.1'))

    f_o = open(fname_o,'r')
    f_c = open(fname_c,'r')

    fout_o = fileName + '_open_v'
    fout_c = fileName + '_close_v'
    f1 = open(fout_o, 'w')
    f2 = open(fout_c, 'w')
    
    customTable(f_o, f1)
    customTable(f_c, f2)

    # FOR QA PURPOSES
    tp.report(FILE, par=('disp','wall_thick'),
	      results=('react.1'), file='elbowtest.psr',
	      variations=OFF)