1 // 1) compile with "g++ pend.cpp -o pend.exe"
2 // 2) launch "gmsh pend.exe"
3 
4 #include <math.h>
5 #include <stdio.h>
6 #include "onelab.h"
7 
exportMsh(const std::string & path,double le1,double le2)8 void exportMsh(const std::string &path, double le1, double le2)
9 {
10   FILE *mshFile = fopen((path + "pend.msh").c_str(),"w");
11   if(!mshFile) return;
12   fprintf(mshFile, "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n");
13   fprintf(mshFile, "$Nodes\n3\n1 0 0 0\n2 0 %f 0\n3 0 %f 0\n$EndNodes\n",
14           -le1, -le1-le2);
15   fprintf(mshFile, "$Elements\n3\n1 1 2 0 1 1 2\n2 1 2 0 1 2 3\n3 15 2 0 2 3\n"
16           "$EndElements\n");
17   fclose(mshFile);
18 }
19 
exportMshOpt(const std::string & path)20 void exportMshOpt(const std::string &path)
21 {
22   FILE *optFile = fopen((path + "pend.msh.opt").c_str(), "w");
23   if(!optFile) return;
24   fprintf(optFile, "n = PostProcessing.NbViews - 1;\n");
25   fprintf(optFile, "If(n >= 0)\nView[n].ShowScale = 0;\nView[n].VectorType = 5;\n");
26   fprintf(optFile, "View[n].ExternalView = 0;\nView[n].DisplacementFactor = 1 ;\n");
27   fprintf(optFile, "View[n].PointType = 1;\nView[n].PointSize = 5;\n");
28   fprintf(optFile, "View[n].LineWidth = 2;\nEndIf\n");
29   fclose(optFile);
30 }
31 
exportIter(const std::string & path,int iter,double t,double x1,double y1,double x2,double y2)32 void exportIter(const std::string &path, int iter, double t, double x1, double y1,
33                 double x2, double y2)
34 {
35   FILE *mshFile = fopen((path + "pend.msh").c_str(), "a");
36   if(!mshFile) return;
37   fprintf(mshFile, "$NodeData\n1\n\"motion\"\n1\n\t%f\n3\n\t%d\n3\n", t, iter);
38   fprintf(mshFile, "\t3\n\t1 0 0 0\n\t2 %f %f 0\n\t3 %f %f 0\n$EndNodeData\n",
39           x1, y1, x2, y2);
40   fclose(mshFile);
41 }
42 
defineNumber(onelab::client * c,const std::string & name,double value,const std::map<std::string,std::string> & attributes)43 double defineNumber(onelab::client *c, const std::string &name, double value,
44                     const std::map<std::string, std::string> &attributes)
45 {
46   std::vector<onelab::number> ns;
47   c->get(ns, name);
48   if(ns.empty()){ // define new parameter
49     onelab::number n(name, value);
50     if(attributes.size()) n.setAttributes(attributes);
51     c->set(n);
52     return value;
53   }
54   // return value from server
55   return ns[0].getValue();
56 }
57 
setNumber(onelab::client * c,const std::string & name,double value,double min=0,double max=0,bool visible=true)58 void setNumber(onelab::client *c, const std::string &name, double value,
59                double min=0, double max=0, bool visible=true)
60 {
61   onelab::number n(name, value);
62   n.setMin(min);
63   n.setMax(max);
64   n.setVisible(visible);
65   c->set(n);
66 }
67 
addNumberChoice(onelab::client * c,const std::string & name,double choice)68 void addNumberChoice(onelab::client *c, const std::string &name, double choice)
69 {
70   std::vector<onelab::number> ns;
71   c->get(ns, name);
72   if(ns.size()){
73     std::vector<double> choices = ns[0].getChoices();
74     choices.push_back(choice);
75     ns[0].setChoices(choices);
76     c->set(ns[0]);
77   }
78 }
79 
main(int argc,char ** argv)80 int main(int argc, char **argv)
81 {
82   std::string name, address;
83   for(int i = 0; i < argc; i++){
84     if(std::string(argv[i]) == "-onelab" && i + 2 < argc){
85       name = std::string(argv[i + 1]);
86       address = std::string(argv[i + 2]);
87     }
88   }
89 
90   if(name.empty() || address.empty()) return 1;
91 
92   onelab::remoteNetworkClient *c = new onelab::remoteNetworkClient(name, address);
93 
94   std::string action;
95   std::vector<onelab::string> ns;
96   c->get(ns, name + "/Action");
97   if(ns.size()) action = ns[0].getValue();
98 
99   // prevent automatic Gmsh model reload & meshing
100   onelab::number n("IsMetamodel", 1);
101   c->set(n);
102 
103   std::string path(argv[0]);
104   int islash = (int)path.find_last_of("/\\");
105   if(islash > 0)
106     path = path.substr(0, islash + 1);
107   else
108     path = "";
109 
110   double g = 9.8; // acceleration of gravity
111   double m = 0.3; // mass of pendulum balls
112 
113   std::map<std::string, std::string> attr;
114 
115   double l = defineNumber(c, "Geom/arm length [m]", 1.0, attr);
116   double time = defineNumber(c, "Dyna/time [s]", 0., attr);
117   double dt = defineNumber(c, "Dyna/time step [s]", 0.001, attr);
118   double tmax = defineNumber(c, "Dyna/max time [s]", 20, attr);
119   double refresh = defineNumber(c, "Dyna/refresh interval [s]", 0.05, attr);
120   attr["Highlight"] = "Pink";
121   double theta0 = defineNumber(c, "Init/initial theta angle [deg]", 10, attr);
122   double phi0 = defineNumber(c, "Init/initial phi angle [deg]", 180, attr);
123 
124   // we're done if we are not in the compute phase
125   if(action != "compute"){
126     delete c;
127     return 0;
128   }
129 
130   double l1 = l;
131   double l2 = l;
132   double m1 = m;
133   double m2 = m;
134   double theta = theta0 / 180.*M_PI;
135   double phi = phi0 / 180.*M_PI;
136   double theta_dot = 0.0;
137   double phi_dot = 0.0;
138   double refr = 0.0;
139   int iter = 0;
140   time = 0.0;
141 
142   while (time < tmax){
143     double delta = phi - theta;
144     double sdelta = sin(delta);
145     double cdelta = cos(delta);
146     double theta_dot_dot = ( m2*l1*(theta_dot*theta_dot)*sdelta*cdelta
147                              + m2*g*sin(phi)*cdelta
148                              + m2*l2*(phi_dot*phi_dot)*sdelta
149                              - (m1+m2)*g*sin(theta) );
150     theta_dot_dot /= ( (m1+m2)*l1 - m2*l1*(cdelta*cdelta) );
151 
152     double phi_dot_dot = ( -m2*l2*(phi_dot*phi_dot)*sdelta*cdelta
153                            + (m1+m2)*(g*sin(theta)*cdelta
154                                       - l1*(theta_dot*theta_dot)*sdelta
155                                       - g*sin(phi)) );
156     phi_dot_dot /= ( (m1+m2)*l2 - m2*l2*(cdelta*cdelta) );
157 
158     theta_dot += theta_dot_dot*dt;
159     phi_dot += phi_dot_dot*dt;
160     theta += theta_dot*dt;
161     phi += phi_dot*dt;
162 
163     double x1 =  l1*sin(theta);
164     double y1 = -l1*cos(theta);
165     double x2 =  l1*sin(theta) + l2*sin(phi);
166     double y2 = -l1*cos(theta) - l2*cos(phi);
167 
168     time += dt;
169     refr += dt;
170 
171     exportMshOpt(path);
172 
173     if(refr >= refresh){
174       refr = 0;
175       setNumber(c, name + "/Progress", time, 0, tmax, false);
176       setNumber(c, "Dyna/time [s]", time);
177       setNumber(c, "Solu/phi", phi);
178       addNumberChoice(c, "Solu/phi", phi);
179       setNumber(c, "Solu/theta",  theta);
180       addNumberChoice(c, "Solu/theta", theta);
181       setNumber(c, "Solu/phi dot", phi_dot);
182       addNumberChoice(c, "Solu/phi dot", phi_dot);
183       setNumber(c, "Solu/theta dot", theta_dot);
184       addNumberChoice(c, "Solu/theta dot", theta_dot);
185 
186       // ask Gmsh to refresh
187       onelab::string s("Gmsh/Action", "refresh");
188       c->set(s);
189 
190       // stop if we are asked to (by Gmsh)
191       c->get(ns, name + "/Action");
192       if(ns.size() && ns[0].getValue() == "stop") break;
193 
194       exportMsh(path, l1, l2);
195       exportIter(path, iter, time, x1, y1+l1, x2, y2+l1+l2);
196       c->sendMergeFileRequest(path + "pend.msh");
197       iter += 1;
198     }
199   }
200 
201   setNumber(c, name + "/Progress", 0, 0, tmax, false);
202 
203   delete c;
204 
205   return 0;
206 }
207