added skeleton of FOLD input

ex/simulate_FOLD.py

+#!/usr/bin/env python
+from __future__ import division
+from numpy import *
+import argparse
+from pyframe3dd.frame3dd import write_frame3dd_file, read_lowest_mode, read_frame3dd_displacements, compute_mass
+from pyframe3dd.util import magnitudes, close
+import subprocess
+
+def run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads):
+	write_frame3dd_file(nodes,global_args,beam_sets,constraints,loads)
+	cmd = ["frame3dd", "-i",global_args['frame3dd_filename']+'.csv']
+	if args.quiet: cmd.append("-q")
+	print ' '.join(cmd)
+	subprocess.call(cmd)
+
+def clean_up_frame3dd(filename):
+	#Delete files generated by frame3dd
+	files = [filename+end for end in ["_out.csv",".csv.out",".csv.plt",".csv.if01",".csv"]]
+	subprocess.call(["rm"]+files)
+
+def build(args):
+	#return nodes,rods as numpy arrays
+	#open fold file based on base_filename argument
+	#parse into dictionary
+	#get vertices from "vertices_coords"
+	#get edges from "edges_vertices"
+	#return numpy arrays
+
+	#once the above is happy, you can check if faces_vertices exists, and if so, get edges from faces_vertices by removing duplicates.
+	#once that is happy, we can code beam templates from the tomohiro tachi paper to apply to each face.
+	return args.l*array([[0,0,0],[0,0,1.]]), array([[0,1]])
+
+def run_simulation(args):
+	#set up simulation
+	nodes,rods = build(args)
+	global_args = {
+		'n_modes':args.n_modes,'length_scaling':args.length_scaling,
+		'node_radius':zeros(shape(nodes)[0]),'frame3dd_filename':args.base_filename+"_frame3dd"
+	}
+	clean_up_frame3dd(global_args['frame3dd_filename'])
+	beam_sets = [
+	(rods,{'E':args.E,'nu':args.nu,'rho':args.rho,'cross_section':'circular','d1':args.d,'th':args.t,'roll':0.,'loads':[],'beam_divisions':args.bd,'prestresses':[]})
+	]
+	
+	constraints =      [{'node':0,'DOF':dof,'value':0} for dof in [0,1,2,5]]
+	constraints.extend([{'node':1,'DOF':dof,'value':0} for dof in [0,1]])
+
+	loads = [{'node':1,'DOF':2,'value':-args.force}]
+
+	run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
+	results = {}
+	results['beam_mass'] = compute_mass(nodes,beam_sets)
+	results['fundamental_frequency'] = read_lowest_mode(global_args['frame3dd_filename']+'.csv')
+	Ro = .5*args.d; t = args.t; a = pi*(Ro**2 - (Ro-t)**2)
+	results['stress'] = args.force/a
+	return results
+
+def find_stability_threshold(args):
+	#out loop of simulations to determine the buckling load
+	lower = 0 #lower bound
+	upper = 10*args.force_res #initial upper bound before bracketing
+	bracketed=False
+	#actually not necessary, but fun to have the unloaded frequency
+	args.force = lower
+	res = run_simulation(args)	
+	freqs = [res['fundamental_frequency']]
+	forces = [args.force]
+
+	i = 0
+	while not bracketed:
+		print lower,upper,bracketed,res['fundamental_frequency']
+		args.force = upper
+		res = run_simulation(args); i += 1
+		if res['fundamental_frequency']<0:
+			bracketed=True
+		else:
+			freqs.append(res['fundamental_frequency'])
+			forces.append(args.force)
+			lower = upper
+			upper = 2*upper
+	while (upper-lower > args.force_res):
+		print lower,upper,bracketed
+		args.force = .5*(upper+lower)
+		res = run_simulation(args); i += 1
+		if res['fundamental_frequency']>0:
+			freqs.append(res['fundamental_frequency'])
+			forces.append(args.force)
+			lower = .5*(upper+lower)
+		else:
+			upper = .5*(upper+lower)
+	return forces,freqs,res
+
+
+if __name__ == '__main__':
+	parser = argparse.ArgumentParser()
+	parser.add_argument('-M','--mode',choices=('simulate','search'), required=True)
+	parser.add_argument('-Q','--quiet',action='store_true',help='Whether to suppress frame3dd output')
+	parser.add_argument("-f","--force",  type=double, default=10e90, help="force to apply")
+	parser.add_argument("-fr","--force_res", type=double, default=.01, help="Final resolution of force for search mode")
+
+	parser.add_argument("-d","--d", type=double, default=.01, help="diameter of tube")
+	parser.add_argument("-t","--t", type=double, default=.001, help="wall thickness of tube")
+	parser.add_argument("-l","--l", type=double, default=.5, help="length of tube")
+
+	parser.add_argument("-bd","--bd", type=int, default=1, help='how many divisions for each rod, useful in buckling analysis')
+	parser.add_argument("-E","--E", type=double, default=10e9, help="Young's Modulus of laminate")
+	parser.add_argument("-nu","--nu", type=double, default=.35, help="Poisson Ratio")
+	parser.add_argument("-base_filename","--base_filename", default='buckle', help="Base filename for segments and frame3dd")
+	parser.add_argument("-rho","--rho",type=double,default=1600.,help='density of beam material, kg/m^3')
+	parser.add_argument("-n_modes","--n_modes",type=int,default=4,help='number of modes to compute')
+	parser.add_argument("-ls","--length_scaling", type=double, default=1., help="Scale factor to keep numbers commesurate")
+	args = parser.parse_args()
+
+	if args.mode=='search':
+		forces,freqs,last_res = find_stability_threshold(args)
+		print "Fundamental frequency: %.3f Hz"%(freqs[-1])
+		print "Critical force: %.3f N"%(forces[-1])
+		print "Critical stress: %.3f MPa"%(last_res['stress']/1e6)
+	elif args.mode=='simulate':
+		res = run_simulation(args)
+		print "Fundamental frequency: %.3f Hz"%res['fundamental_frequency']
+		print "Stress: %.3f MPa"%(res['stress']/1e6)
+	else:
+		assert(0) #should not be here
+
+
+