1#!/usr/bin/env python 2#coding=utf-8 3 4# 1) launch "gmsh pend.py" 5# 2) there is no 2... :-) 6 7import onelab 8import math, os 9 10c = onelab.client(__file__) 11 12def exportMsh(le1,le2): 13 mshFile = open(c.getPath("pend.msh"), 'w') 14 mshFile.write('$MeshFormat\n2.2 0 8\n$EndMeshFormat\n') 15 mshFile.write('$Nodes\n3\n1 0 0 0\n2 0 %s 0\n3 0 %s 0\n$EndNodes\n' %(-le1, -le1-le2)) 16 mshFile.write('$Elements\n3\n1 1 2 0 1 1 2\n2 1 2 0 1 2 3\n3 15 2 0 2 3\n$EndElements\n') 17 mshFile.close() 18 19def exportMshOpt(): 20 optFile = open(c.getPath("pend.msh.opt"),'w') 21 optFile.write('n = PostProcessing.NbViews - 1;\n') 22 optFile.write('If(n >= 0)\nView[n].ShowScale = 0;\nView[n].VectorType = 5;\n') 23 optFile.write('View[n].ExternalView = 0;\nView[n].DisplacementFactor = 1 ;\n') 24 optFile.write('View[n].PointType = 1;\nView[n].PointSize = 5;\n') 25 optFile.write('View[n].LineWidth = 2;\nEndIf\n') 26 optFile.close() 27 28def exportIter(iter,t,x1,y1,x2,y2): 29 mshFile = open(c.getPath("pend.msh"),'a') 30 mshFile.write('$NodeData\n1\n"motion"\n1\n\t%f\n3\n\t%d\n3\n' % (t, iter)) 31 mshFile.write('\t3\n\t1 0 0 0\n\t2 %f %f 0\n\t3 %f %f 0\n$EndNodeData\n' %(x1,y1,x2,y2)) 32 mshFile.close() 33 34g = 9.8 # acceleration of gravity 35m = 0.3 # mass of pendulum balls 36 37l = c.defineNumber('Geom/arm length [m]', value=1.0) 38time = c.defineNumber('Dyna/time [s]', value=0.0) 39dt = c.defineNumber('Dyna/time step [s]', value=0.001) 40tmax = c.defineNumber('Dyna/max time [s]', value=20) 41refresh = c.defineNumber('Dyna/refresh interval [s]', value=0.1) 42theta0 = c.defineNumber('Init/initial theta angle [deg]', value=10, 43 attributes={'Highlight':'Pink'}) 44phi0 = c.defineNumber('Init/initial phi angle [deg]', value=180, 45 attributes={'Highlight':'Pink'}) 46 47# we're done if we are in the "check" phase 48if c.action == 'check' : 49 exit(0) 50 51l1 = l; 52l2 = l; 53m1 = m; 54m2 = m; 55theta = theta0 / 180.*math.pi; 56phi = phi0 / 180.*math.pi; 57theta_dot = 0.0 58phi_dot = 0.0 59refr = 0.0 60iter = 0 61time = 0.0 62 63while (time < tmax): 64 delta = phi - theta 65 sdelta = math.sin(delta) 66 cdelta = math.cos(delta) 67 theta_dot_dot = ( m2*l1*(theta_dot**2.0)*sdelta*cdelta 68 + m2*g*math.sin(phi)*cdelta 69 + m2*l2*(phi_dot**2.0)*sdelta 70 - (m1+m2)*g*math.sin(theta) ) 71 theta_dot_dot /= ( (m1+m2)*l1 - m2*l1*(cdelta)**2.0 ) 72 73 phi_dot_dot = ( -m2*l2*(phi_dot**2.0)*sdelta*cdelta 74 + (m1+m2)*(g*math.sin(theta)*cdelta 75 - l1*(theta_dot**2.0)*sdelta 76 - g*math.sin(phi)) ) 77 phi_dot_dot /= ( (m1+m2)*l2 - m2*l2*(cdelta)**2.0 ) 78 79 theta_dot = theta_dot + theta_dot_dot*dt 80 phi_dot = phi_dot + phi_dot_dot*dt 81 82 theta = theta + theta_dot*dt 83 phi = phi + phi_dot*dt 84 85 x1 = l1*math.sin(theta) 86 y1 = -l1*math.cos(theta) 87 x2 = l1*math.sin(theta) + l2*math.sin(phi) 88 y2 = -l1*math.cos(theta) - l2*math.cos(phi) 89 90 time += dt 91 refr += dt 92 93 exportMshOpt() 94 95 if refr >= refresh: 96 refr = 0 97 c.setNumber(c.name + '/Progress', value=time, min=0, max=tmax, visible=0) 98 c.setNumber('Dyna/time [s]', value=time) 99 c.setNumber('Solu/phi', value=phi) 100 c.addNumberChoice('Solu/phi', phi) 101 c.setNumber('Solu/theta', value=theta) 102 c.addNumberChoice('Solu/theta', theta) 103 c.setNumber('Solu/phi dot', value=phi_dot) 104 c.addNumberChoice('Solu/phi dot', phi_dot) 105 c.setNumber('Solu/theta dot', value=theta_dot) 106 c.addNumberChoice('Solu/theta dot', theta_dot) 107 108 # ask Gmsh to refresh 109 c.setString('Gmsh/Action', value='refresh') 110 111 # stop if we are asked to (by Gmsh) 112 if(c.getString(c.name + '/Action') == 'stop'): 113 break; 114 115 exportMsh(l1, l2) 116 exportIter(iter, time, x1, y1+l1, x2, y2+l1+l2) 117 c.mergeFile(c.checkPath('pend.msh')) 118 iter += 1 119 120c.setNumber(c.name + '/Progress', value=0) 121