simulate_FOLD.py 8.87 KB
Newer Older
Sam Calisch's avatar
Sam Calisch committed
1 2 3 4 5 6 7
#!/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
Soma Mitra-Behura's avatar
Soma Mitra-Behura committed
8 9
import json
import sys
Sam Calisch's avatar
Sam Calisch committed
10

11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
def plot_connections(args,nodes,connect,constraints=None,loads=None):
    #for debug only, this is slow!
    import matplotlib as mpl
    from mpl_toolkits.mplot3d import Axes3D
    import numpy as np
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    with open(args.base_filename+".fold") as data_file:
        fold = json.load(data_file)
    verts = args.l*asarray(fold["vertices_coords"])
    ax.plot(nodes[:,0],nodes[:,1],nodes[:,2],ls="",marker="o")
    #print loads
    #print len(loads)
    if constraints==None:
    	pass
    elif len(constraints)>0:
    	for cons in constraints:
    		node=cons['node']
    		dof=cons['DOF']
    		val=cons['value']
    		# ax.text(verts[cons['node']][0],verts[cons['node']][1],verts[cons['node']][2],cons['DOF'])
    		#print dof
    		if dof==0:
    			ax.quiver((verts[node][0]-0.1),(verts[node][1]),(verts[node][2]),(1),(0),(0),length=0.1,normalize=True)
    		if dof==1:
				ax.quiver((verts[node][0]),(verts[node][1]-0.1),(verts[node][2]),(0),(1),(0),length=0.1,normalize=True)
    		if dof==2:
				ax.quiver((verts[node][0]),(verts[node][1]),(verts[node][2]-0.1),(0),(0),(1),length=0.1,normalize=True)
	#if loads==None:
	#	print 'got here'
    #	pass
    if len(loads)>0:
    	#print 'got here'
    	for load in loads:
    		node=load['node']
    		dof=load['DOF']
    		print dof
    		if dof==0:
    			ax.quiver((verts[node][0]-0.1),(verts[node][1]),(verts[node][2]),(1),(0),(0),length=0.1,normalize=True,color='r')
    		if dof==1:
				ax.quiver((verts[node][0]),(verts[node][1]-0.1),(verts[node][2]),(0),(1),(0),length=0.1,normalize=True,color='r')
    		if dof==2:
				ax.quiver((verts[node][0]),(verts[node][1]),(verts[node][2]-0.1),(0),(0),(1),length=0.1,normalize=True,color='r')
    for seg in connect:
        ax.plot(nodes[seg,0], nodes[seg,1], nodes[seg,2], c='g')
    plt.show()

Sam Calisch's avatar
Sam Calisch committed
59 60 61 62 63 64 65 66 67 68 69 70 71
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):
Soma Mitra-Behura's avatar
Soma Mitra-Behura committed
72
	#return nodes,rods as numpy arrays	#open fold file based on base_filename argument
Sam Calisch's avatar
Sam Calisch committed
73 74 75 76 77 78 79
	#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.
Soma Mitra-Behura's avatar
Soma Mitra-Behura committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
	with open(args.base_filename+".fold") as data_file:
		fold = json.load(data_file)
	verts = asarray(fold["vertices_coords"])
	edges = asarray(fold["edges_vertices"])
	faces = asarray(fold["faces_vertices"])
	if args.facet_model == 'N4B5':
		for face in faces:
			try:
				p4 = face[3]
				edges = vstack([edges, array([face[0],face[2]])])
			except:
				print 'One of input faces is a triangle.'
	if args.facet_model == 'N4B6':
		for face in faces:
			try:
				edges = vstack([edges, [array([face[0],face[2]]),array([face[1],face[3]])]])
			except:
				print 'One of input faces is a triangle.'
	if args.facet_model == 'N5B8':
		for face in faces:
			try:
				p0 = verts[face[0]]
				p1 = verts[face[1]]
				p2 = verts[face[2]]
				p3 = verts[face[3]]
				xs = [p0[0],p1[0],p2[0],p3[0]]
				ys = [p0[1],p1[1],p2[1],p3[1]]
				zs = [p0[2],p1[2],p2[2],p3[2]]
				verts = vstack([verts,array([mean(xs),mean(ys),mean(zs)])])
				edges = vstack([edges,[array([face[0],len(verts)-1]),array([face[1],len(verts)-1]),array([face[2],len(verts)-1]),array([face[3],len(verts)-1])]])
			except:
				print 'One of input faces is a triangle.'
	return args.l*verts, edges
	# return args.l*array([[0,0,0],[0,0,1.]]), array([[0,1]])
Sam Calisch's avatar
Sam Calisch committed
114

115 116 117 118 119 120 121 122 123 124 125 126
def build_constraints():
	constraints =      [{'node':0,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5]]
	constraints.extend([{'node':1,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5]])
	constraints.extend([{'node':3,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5]])
	constraints.extend([{'node':4,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5]])
	constraints.extend([{'node':5,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5]])
	return constraints

def build_loads():
	loads = [{'node':2,'DOF':0,'value':-args.force}]
	return loads

Sam Calisch's avatar
Sam Calisch committed
127 128 129 130 131 132 133 134 135 136 137 138
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':[]})
	]
	
139 140
	constraints = build_constraints()
	print 'These are constraints:', constraints
Sam Calisch's avatar
Sam Calisch committed
141

142 143
	loads = build_loads()
	print 'The are loads', loads
Sam Calisch's avatar
Sam Calisch committed
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190

	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()
191
	parser.add_argument('-M','--mode',choices=('simulate','search','visualize'), required=True)
Sam Calisch's avatar
Sam Calisch committed
192
	parser.add_argument('-Q','--quiet',action='store_true',help='Whether to suppress frame3dd output')
193
	parser.add_argument("-f","--force",  type=double, default=1, help="force to apply")
Sam Calisch's avatar
Sam Calisch committed
194 195 196 197
	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")
Soma Mitra-Behura's avatar
Soma Mitra-Behura committed
198
	parser.add_argument("-l","--l", type=double, default=0.5, help="length of tube") #used to be type=double, default=0.5
Sam Calisch's avatar
Sam Calisch committed
199 200 201 202 203 204 205 206

	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")
Soma Mitra-Behura's avatar
Soma Mitra-Behura committed
207
	parser.add_argument("-fm","--facet_model",choices=('N4B4','N4B5','N4B6','N5B8'),default='N4B4')
Sam Calisch's avatar
Sam Calisch committed
208 209
	args = parser.parse_args()

Soma Mitra-Behura's avatar
Soma Mitra-Behura committed
210 211
	print args.base_filename

Sam Calisch's avatar
Sam Calisch committed
212 213 214 215 216 217 218 219 220
	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)
221 222 223
	elif args.mode=='visualize':
		nodes,rods = build(args)
		plot_connections(args,nodes,rods,build_constraints(),build_loads())
Sam Calisch's avatar
Sam Calisch committed
224 225 226 227 228
	else:
		assert(0) #should not be here