diff --git a/ex/simulate_FOLD.py b/ex/simulate_FOLD.py index 2be3f1cafba4b5cb497cdd3fbaa089b08bd1ae09..db9afd356a1c8ddefd1fab63afab3f6a16bdc5e6 100644 --- a/ex/simulate_FOLD.py +++ b/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