Commit 57b2558a authored by Soma Mitra-Behura's avatar Soma Mitra-Behura

reupdated simulate

parent 5cce00fe
......@@ -5,6 +5,8 @@ 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
import json
import sys
def run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads):
write_frame3dd_file(nodes,global_args,beam_sets,constraints,loads)
......@@ -19,8 +21,7 @@ def clean_up_frame3dd(filename):
subprocess.call(["rm"]+files)
def build(args):
#return nodes,rods as numpy arrays
#open fold file based on base_filename argument
#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"
......@@ -28,7 +29,61 @@ def build(args):
#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]])
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':
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)])])
# 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 run_simulation(args):
#set up simulation
......@@ -42,10 +97,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,5]]
constraints.extend([{'node':1,'DOF':dof,'value':0} for dof in [0,1]])
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]])
loads = [{'node':1,'DOF':2,'value':-args.force}]
loads = [{'node':3,'DOF':2,'value':-args.force}]
run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
results = {}
......@@ -100,7 +156,7 @@ if __name__ == '__main__':
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("-l","--l", type=double, default=0.5, help="length of tube") #used to be type=double, default=0.5
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")
......@@ -109,8 +165,11 @@ if __name__ == '__main__':
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")
parser.add_argument("-fm","--facet_model",choices=('N4B4','N4B5','N4B6','N5B8'),default='N4B4')
args = parser.parse_args()
print args.base_filename
if args.mode=='search':
forces,freqs,last_res = find_stability_threshold(args)
print "Fundamental frequency: %.3f Hz"%(freqs[-1])
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment