updated simulate fold python file

ex/simulate_FOLD.py

@@ -8,6 +8,54 @@ import subprocess
 import json
 import sys
 
+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()
+
 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']
@@ -32,31 +80,19 @@ def build(args):
 	with open(args.base_filename+".fold") as data_file:
 		fold = json.load(data_file)
 	verts = asarray(fold["vertices_coords"])
-	print verts
 	edges = asarray(fold["edges_vertices"])
-	print edges
 	faces = asarray(fold["faces_vertices"])
-	print faces
-	print len(faces)
 	if args.facet_model == 'N4B5':
 		for face in faces:
 			try:
-				print 'Current face is', face, 'and type is', type(face)
 				p4 = face[3]
-				print 'New edge is', array([face[0],face[2]])
-				# edges = append(edges, array([face[0],face[2]]))
-				# concatenate(edges, array([face[0],face[2]]))
 				edges = vstack([edges, array([face[0],face[2]])])
-				# edges.append(array([face[0],face[2]]))
-				print 'Updated edges list is', edges
-				sys.exit(0)
 			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]])]])
-				print edges
 			except:
 				print 'One of input faces is a triangle.'
 	if args.facet_model == 'N5B8':
@@ -70,21 +106,24 @@ def build(args):
 				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)])])
-				# verts.append([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])]])
-				print verts
-				print edges
-				# edges.append([face[0],len(verts)-1])
-				# edges.append([face[1],len(verts)-1])
-				# edges.append([face[2],len(verts)-1])
-				# edges.append([face[3],len(verts)-1])
 			except:
 				print 'One of input faces is a triangle.'
-	# print verts, edges
-	# sys.exit(0)
 	return args.l*verts, edges
 	# return args.l*array([[0,0,0],[0,0,1.]]), array([[0,1]])
 
+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
+
 def run_simulation(args):
 	#set up simulation
 	nodes,rods = build(args)
@@ -97,11 +136,11 @@ def run_simulation(args):
 	(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]]
-	constraints.extend([{'node':1,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5]])
-	constraints.extend([{'node':2,'DOF':dof,'value':0} for dof in [0,1,2]])
+	constraints = build_constraints()
+	print 'These are constraints:', constraints
 
-	loads = [{'node':3,'DOF':2,'value':-args.force}]
+	loads = build_loads()
+	print 'The are loads', loads
 
 	run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
 	results = {}
@@ -149,9 +188,9 @@ def find_stability_threshold(args):
 
 if __name__ == '__main__':
 	parser = argparse.ArgumentParser()
-	parser.add_argument('-M','--mode',choices=('simulate','search'), required=True)
+	parser.add_argument('-M','--mode',choices=('simulate','search','visualize'), 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("-f","--force",  type=double, default=1, 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")
@@ -179,6 +218,9 @@ if __name__ == '__main__':
 		res = run_simulation(args)
 		print "Fundamental frequency: %.3f Hz"%res['fundamental_frequency']
 		print "Stress: %.3f MPa"%(res['stress']/1e6)
+	elif args.mode=='visualize':
+		nodes,rods = build(args)
+		plot_connections(args,nodes,rods,build_constraints(),build_loads())
 	else:
 		assert(0) #should not be here