diff --git a/flexure/CMD b/flexure/CMD
new file mode 100644
index 0000000000000000000000000000000000000000..25049fe8183db196e91084e01ee931d177bfe89d
--- /dev/null
+++ b/flexure/CMD
@@ -0,0 +1 @@
+python flexure_stiffness.py -Q -M simulate -f 15 -flexure_type cyclic -w .0007 -t .0015 -l .010
diff --git a/flexure/flexure_stiffness.py b/flexure/flexure_stiffness.py
new file mode 100644
index 0000000000000000000000000000000000000000..ba1910f406dad9a26dfe4470c2b8366878873b4c
--- /dev/null
+++ b/flexure/flexure_stiffness.py
@@ -0,0 +1,183 @@
+#!/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 plot_connections(nodes,beamsets):
+ #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()
+ cmap = plt.cm.get_cmap('Dark2', len(beamsets)+1)
+ ax = fig.gca(projection='3d')
+ for i,beamset in enumerate(beamsets):
+ for seg in beamset:
+ ax.plot(nodes[seg,0], nodes[seg,1], nodes[seg,2], c=cmap(i))
+ 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']
+ 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
+ dxy = args.attach_radius/sqrt(2)
+
+ if args.flexure_type == 'mirrored':
+ nodes = args.l*array([[0,0,0],[0,1,0],[-1,1,0],[1,0,0],[1,-1,0]]) #one side of flexure plate
+ beams = array([[0,1],[1,2],[0,3],[3,4]])
+ nodes += array([dxy,dxy,0]) #offset by attachment radius
+
+ elif args.flexure_type == 'cyclic':
+ nodes = args.l*array([[0,0,0],[0,1,0],[-1,1,0]]) + array([dxy,dxy,0])
+ beams = array([[0,1],[1,2],[3,4]])
+ nodes = vstack((nodes, args.l*array([[1,1,0],[1,0,0]]) + array([dxy,-dxy,0])))
+
+ nodes = vstack((nodes,array([[-n[0],-n[1],0] for n in nodes]))) #append reflection
+ beams = vstack((beams, beams + 5))
+ z_os = array([0,0,.5*args.sep])
+ nodes = vstack((nodes + z_os, nodes - z_os))
+ beams = vstack((beams, beams + 10))
+
+ solid_nodes = array([ [dxy,-dxy,.5*args.sep],[-dxy,dxy,.5*args.sep],[dxy,-dxy,-.5*args.sep],[-dxy,dxy,-.5*args.sep] ])
+ nodes = vstack((nodes, solid_nodes))
+ solid_beams = array([
+ [0,5],[0,10],[5,15],[10,15],#[0,15],[5,10],
+ [0,20],[0,21],[5,20],[5,21],[20,21],
+ [10,22],[10,23],[15,22],[15,23],[22,23],
+ [20,22],[21,23],#[20,23],[21,22],
+ [0,22],[0,23],[5,22],[5,23],
+ [10,20],[10,21],[15,20],[15,21]
+ ])
+ if args.flexure_type == 'cyclic':
+ beams = vstack((beams, array([[4,20],[9,21],[14,22],[19,23]])))
+ return nodes, beams, solid_beams
+
+def run_simulation(args):
+ #set up simulation
+ nodes,beams,solid_beams = build(args)
+ global_args = {
+ 'n_modes':args.n_modes,'length_scaling':args.length_scaling,'exagerration':10,
+ 'node_radius':zeros(shape(nodes)[0]),'frame3dd_filename':args.base_filename+"_frame3dd"
+ }
+ clean_up_frame3dd(global_args['frame3dd_filename'])
+ beam_sets = [
+ (beams,{'E':args.E,'nu':args.nu,'rho':args.rho,'cross_section':'rectangular','d2':args.w,'d1':args.t,'roll':0.,'loads':[],'beam_divisions':args.bd,'prestresses':[]}),
+ (solid_beams,{'E':10*args.E,'nu':args.nu,'rho':args.rho,'cross_section':'rectangular','d1':args.w,'d2':args.t,'roll':0.,'loads':[],'beam_divisions':args.bd,'prestresses':[]})
+ ]
+ if args.flexure_type == 'mirrored':
+ fixed_nodes = [2,4,7,9,12,14,17,19]
+ elif args.flexure_type == 'cyclic':
+ fixed_nodes = [2,3,7,8,12,13,17,18]
+ constraints = [{'node':node,'DOF':dof,'value':0} for dof in [0,1,2,3,4,5] for node in fixed_nodes]
+
+
+ for force_dof in [0,1,2]:
+ loaded_nodes = [0,5,10,15]
+ loads = [{'node':n,'DOF':force_dof,'value':args.force/len(loaded_nodes)} for n in loaded_nodes]
+
+ run_frame3dd(args,nodes,global_args,beam_sets,constraints,loads)
+ results = {}
+ results['beam_mass'] = compute_mass(nodes,beam_sets)
+ disps = read_frame3dd_displacements(global_args['frame3dd_filename'])
+ force_disp = average(disps[loaded_nodes,force_dof])
+ print "Degree of freedom: %d"%force_dof
+ print "Force applied: %.1f N"%(args.force)
+ print "Calculated displacement: %.1f microns"%(force_disp*1e6)
+
+ #todo: read average displacement of loaded nodes
+ #todo: plot displacements vs. design parameters
+
+
+ #results['fundamental_frequency'] = read_lowest_mode(global_args['frame3dd_filename']+'.csv')
+ 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', 'visualize'), required=True)
+ parser.add_argument('-flexure_type','--flexure_type',choices=('cyclic','mirrored'), required=True)
+ parser.add_argument('-Q','--quiet',action='store_true',help='Whether to suppress frame3dd output')
+ 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("-w","--w", type=double, default=.0005, help="width of flexure (m)")
+ parser.add_argument("-t","--t", type=double, default=.0023, help="thickness of flexure material (m)")
+ parser.add_argument("-l","--l", type=double, default=.0068, help="length of flexure segment (m)")
+ parser.add_argument("-attach_radius","--attach_radius", type=double, default=.0043, help="distance from z axis to flexure attachment (m)")
+ parser.add_argument("-sep","--sep", type=double, default=.025, help="flexure plate z separation (m)")
+
+ 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=70e9, help="Young's Modulus of laminate")
+ parser.add_argument("-nu","--nu", type=double, default=.33, 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=2700.,help='density of beam material, kg/m^3')
+ parser.add_argument("-n_modes","--n_modes",type=int,default=0,help='number of dynamic 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)
+ elif args.mode=='visualize':
+ nodes,rods,solid_beams = build(args)
+ plot_connections(nodes,[rods,solid_beams])
+ else:
+ assert(0) #should not be here
+
+
+