StiffnessIO.cpp
Go to the documentation of this file.
1 #include "choreo_task_sequence_planner/utils/StiffnessIO.h"
2 
3 /*
4 * OUTPUT_PATH
5 * return path for output files using specific output directory
6 */
7 void StiffnessIO::OutputPath(const char *fname, char fullpath[], const int len, char *default_outdir, int verbose)
8 {
9  int res;
10  assert(fname != NULL);
11 
12  const char *outdir;
13  outdir = getenv("FRAMEFAB_OUTDIR");
14  if (outdir == NULL)
15  {
16  if (default_outdir == NULL)
17  outdir = TempDir();
18  else
19  outdir = default_outdir;
20  }
21 
22  res = sprintf(fullpath, "%s%c%s", outdir, '\\', fname);
23 
24  if (res > len)
25  {
26  printf("ERROR: unable to construct output filename: overflow.\n");
27  exit(16);
28  }
29 
30  if (verbose)
31  {
32  printf("StiffnessIO->OuputPath: Output file path generated: %s\n", fullpath); /* debug */
33  }
34 }
35 
36 
37 /*
38 * GnuPltStaticMesh - create mesh data of deformed and undeformed mesh 25/Nov/2015
39 * use gnuplot
40 * useful gnuplot options: unset xtics ytics ztics border view key
41 * This function illustrates how to read the internal force output data file.
42 * The internal force output data file contains all the information required
43 * to plot deformed meshes, internal axial force, internal shear force, internal
44 * torsion, and internal bending moment diagrams.
45 */
46 void StiffnessIO::GnuPltStaticMesh(
47  const char *fpath,
48  const char *meshpath, const char *plotpath,
49  VX &D,
50  double exagg_static, float scale,
51  DualGraph *ptr_dualgraph, WireFrame *ptr_frame
52  )
53 {
54  int nN = ptr_dualgraph->SizeOfFaceList();
55  int nE = ptr_dualgraph->SizeOfVertList();
56  WireFrame *ptr_wf = ptr_dualgraph->ptr_frame_;
57  std::vector<WF_edge*> wf_edge_list = *ptr_wf->GetEdgeList();
58  std::vector<DualFace*> dual_face_list = *ptr_dualgraph->GetFaceList();
59 
60  int anlyz = 1;
61 
62  FILE *fpif = NULL, *fpm = NULL;
63  double mx, my, mz; /* coordinates of the frame element number labels */
64  char fnif[FILENMAX], meshfl[FILENMAX],
65  D3 = ' ',
66  errMsg[MAXL],
67  ch = 'a';
68  int sfrv = 0, /* *scanf return value */
69  frel, nx, /* frame element number, number of increments */
70  n1, n2; /* node numbers */
71  float x1, y1, z1, /* coordinates of node n1 */
72  x2, y2, z2; /* coordinates of node n2 */
73  int j = 0, m = 0, n = 0,
74  X = 0, Y = 0, Z = 0,
75  lw = 1; /* line width of deformed mesh */
76  time_t now; /* modern time variable type */
77 
78  (void)time(&now);
79 
80  string title = "Fiber Deform Test";
81  // write gnuplot plotting script commands
82 
83  //for (j = 0; j < nN; j++)
84  //{
85  // // check for three-dimensional frame
86  // if (xyz[j][0] != 0.0) X = 1;
87  // if (xyz[j][1] != 0.0) Y = 1;
88  // if (xyz[j][2] != 0.0) Z = 1;
89  //}
90  //
91  //if ((X && Y && Z) || D3_flag)
92  //{
93  // D3 = ' '; D2 = '#';
94  //}
95  //else
96  //{
97  // D3 = '#'; D2 = ' ';
98  //}
99  // open plotting script file for writing
100  if ((fpm = fopen(plotpath, "w")) == NULL)
101  {
102  sprintf(errMsg, "\n error: cannot open gnuplot script file: %s \n", plotpath);
103  printf(errMsg);
104  exit(23);
105  }
106 
107  // file name for deformed mesh data
108  sprintf(meshfl, "%sf.%03d", meshpath, 0);
109 
110  // write header, plot-setup cmds, node label, and element label data
111 
112  // header & node number & element number labels
113 
114  fprintf(fpm, "# FIBERPRINT FRAME STRUCTUAL ANALYSIS RESULTS GCL@USTC");
115  fprintf(fpm, " VERSION %s \n", std::string(VERSION).c_str());
116  fprintf(fpm, "# %s\n", title.c_str());
117  fprintf(fpm, "# %s", ctime(&now));
118  fprintf(fpm, "# G N U P L O T S C R I P T F I L E \n");
119  /* fprintf(fpm,"# X=%d , Y=%d , Z=%d, D3=%d \n", X,Y,Z,D3_flag); */
120 
121  fprintf(fpm, "set autoscale\n");
122  fprintf(fpm, "unset border\n");
123  fprintf(fpm, "set pointsize 1.0\n");
124  fprintf(fpm, "set xtics; set ytics; set ztics; \n");
125  fprintf(fpm, "unset zeroaxis\n");
126  fprintf(fpm, "unset key\n");
127  fprintf(fpm, "unset label\n");
128  fprintf(fpm, "set size ratio -1 # 1:1 2D axis scaling \n");
129  fprintf(fpm, "# set view equal xyz # 1:1 3D axis scaling \n");
130 
131  fprintf(fpm, "# NODE NUMBER LABELS\n");
132  for (j = 0; j < nN; j++)
133  {
134  int o_id = ptr_dualgraph->v_orig_id(j);
135 
136  fprintf(fpm, "set label ' %d' at %12.4e, %12.4e, %12.4e\n",
137  j + 1, ptr_wf->GetPosition(o_id).x(), ptr_wf->GetPosition(o_id).y(), ptr_wf->GetPosition(o_id).z());
138  }
139 
140  fprintf(fpm, "# ELEMENT NUMBER LABELS\n");
141  for (m = 0; m < nE; m++)
142  {
143  int e_id = ptr_dualgraph->e_orig_id(m);
144  WF_edge *ei = wf_edge_list[e_id];
145  int u = ei->ppair_->pvert_->ID();
146  int v = ei->pvert_->ID();
147  double L = ei->Length();
148 
149  int dual_u = ptr_dualgraph->v_dual_id(u);
150  int dual_v = ptr_dualgraph->v_dual_id(v);
151 
152  mx = 0.5 * (ptr_wf->GetPosition(u).x() + ptr_wf->GetPosition(v).x());
153  my = 0.5 * (ptr_wf->GetPosition(u).y() + ptr_wf->GetPosition(v).y());
154  mz = 0.5 * (ptr_wf->GetPosition(u).z() + ptr_wf->GetPosition(v).z());
155  fprintf(fpm, "set label ' %d' at %12.4e, %12.4e, %12.4e\n",
156  m+1, mx, my, mz);
157  }
158 
159  // 3D plot setup commands
160 
161  fprintf(fpm, "%c set parametric\n", D3);
162  fprintf(fpm, "%c set view 60, 70, %5.2f \n", D3, scale);
163  fprintf(fpm, "%c set view equal xyz # 1:1 3D axis scaling \n", D3);
164  fprintf(fpm, "%c unset key\n", D3);
165  fprintf(fpm, "%c set xlabel 'x'\n", D3);
166  fprintf(fpm, "%c set ylabel 'y'\n", D3);
167  fprintf(fpm, "%c set zlabel 'z'\n", D3);
168  // fprintf(fpm,"%c unset label\n", D3 );
169 
170  // different plot title for each load case
171 
172  fprintf(fpm, "set title \"%s\\n", title.c_str());
173  //fprintf(fpm, "analysis file: %s ", IN_file);
174 
175  if (anlyz)
176  {
177  fprintf(fpm, " deflection exaggeration: %.1f ", exagg_static);
178  }
179 
180  fprintf(fpm, "unset clip; \nset clip one; set clip two\n");
181  fprintf(fpm, "set xyplane 0 \n"); // requires Gnuplot >= 4.6
182 
183  // 2D plot command
184 
185  //fprintf(fpm, "%c plot '%s' u 2:3 t 'undeformed mesh' w lp ",
186  // D2, meshpath);
187  //if (!anlyz) fprintf(fpm, "lw %d lt 1 pt 6 \n", lw);
188  //else fprintf(fpm, "lw 1 lt 5 pt 6, '%s' u 1:2 t 'load case %d of %d' w l lw %d lt 3\n", meshfl, 1, 1, lw);
189 
190  // 3D plot command
191 
192  fprintf(fpm, "%c splot '%s' u 2:3:4 t 'load case %d of %d' w lp ",
193  D3, meshpath, 1, 1);
194  if (!anlyz) fprintf(fpm, " lw %d lt 1 pt 6 \n", lw);
195  else fprintf(fpm, " lw 1 lt 5 pt 6, '%s' u 1:2:3 t 'load case %d of %d' w l lw %d lt 3\n", meshfl, 1, 1, lw);
196 
197  if (anlyz) fprintf(fpm, "pause -1\n");
198 
199  fclose(fpm);
200 
201  // write undeformed mesh data
202 
203  // open the undeformed mesh data file for writing
204  if ((fpm = fopen(meshpath, "w")) == NULL)
205  {
206  sprintf(errMsg, "\n error: cannot open gnuplot undeformed mesh data file: %s\n", meshpath);
207  printf(errMsg);
208  exit(21);
209  }
210 
211  fprintf(fpm, "# FIBERPRINT FRAME STRUCTURAL ANALYSIS RESULTS GCL@USTC");
212  fprintf(fpm, " VERSION %s \n", std::string(VERSION).c_str());
213  fprintf(fpm, "# %s\n", title.c_str());
214  fprintf(fpm, "# %s", ctime(&now));
215  fprintf(fpm, "# U N D E F O R M E D M E S H D A T A (global coordinates)\n");
216  fprintf(fpm, "# Node X Y Z \n");
217 
218  for (m = 0; m < nE; m++)
219  {
220  int e_id = ptr_dualgraph->e_orig_id(m);
221  WF_edge *ei = wf_edge_list[e_id];
222  int u = ei->ppair_->pvert_->ID();
223  int v = ei->pvert_->ID();
224  double L = ei->Length();
225 
226  int dual_u = ptr_dualgraph->v_dual_id(u);
227  int dual_v = ptr_dualgraph->v_dual_id(v);
228 
229  mx = 0.5 * (ptr_wf->GetPosition(u).x() + ptr_wf->GetPosition(v).x());
230  my = 0.5 * (ptr_wf->GetPosition(u).y() + ptr_wf->GetPosition(v).y());
231  mz = 0.5 * (ptr_wf->GetPosition(u).z() + ptr_wf->GetPosition(v).z());
232 
233  fprintf(fpm, "%5d %12.4e %12.4e %12.4e \n",
234  dual_u+1, ptr_wf->GetPosition(u).x(), ptr_wf->GetPosition(u).y(), ptr_wf->GetPosition(u).z());
235  fprintf(fpm, "%5d %12.4e %12.4e %12.4e",
236  dual_v+1, ptr_wf->GetPosition(v).x(), ptr_wf->GetPosition(v).y(), ptr_wf->GetPosition(v).z());
237  fprintf(fpm, "\n\n\n");
238  }
239  fclose(fpm);
240 
241  if (!anlyz) return; // no deformed mesh
242 
243  // write deformed mesh data
244 
245  // open the deformed mesh data file for writing
246  if ((fpm = fopen(meshfl, "w")) == NULL)
247  {
248  sprintf(errMsg, "\n error: cannot open gnuplot deformed mesh data file %s \n", meshfl);
249  printf(errMsg);
250  exit(22);
251  }
252 
253  fprintf(fpm, "# FIBERPRINT FRAME STRUCTURAL ANALYSIS RESULTS GCL@USTC");
254  fprintf(fpm, " VERSION %s \n", std::string(VERSION).c_str());
255  fprintf(fpm, "# %s\n", title.c_str());
256  fprintf(fpm, "# %s", ctime(&now));
257  fprintf(fpm, "# D E F O R M E D M E S H D A T A ");
258  fprintf(fpm, " deflection exaggeration: %.1f\n", exagg_static);
259  fprintf(fpm, "# X-dsp Y-dsp Z-dsp\n");
260 
261 
262  FILE *fpv = fopen("C:/Users/DELL/Desktop/result/vert.txt", "w");
263  FILE *fpl = fopen("C:/Users/DELL/Desktop/result/line.txt", "w");
264 
265  for (m = 0; m < nE; m++)
266  {
267  // write deformed shape data for each element
268 
269  ch = 'a';
270 
271  int e_id = ptr_dualgraph->e_orig_id(m);
272  WF_edge *ei = wf_edge_list[e_id];
273  int u = ei->ppair_->pvert_->ID();
274  int v = ei->pvert_->ID();
275  double L = ei->Length();
276 
277  int dual_u = ptr_dualgraph->v_dual_id(u);
278  int dual_v = ptr_dualgraph->v_dual_id(v);
279 
280  fprintf(fpm, "\n# element %5d \n", m);
281  if (anlyz)
282  {
283  vector<point> beam;
284  GnuPltCubicBentBeam(beam,
285  D, m, ptr_dualgraph, ptr_wf, exagg_static);
286 
287  int nB = beam.size();
288  for (int i = 0; i < nB; ++i)
289  {
290  fprintf(fpm, " %12.4e %12.4e %12.4e\n\n", beam[i].x(), beam[i].y(), beam[i].z());
291  }
292 
293  if (nB != 0)
294  {
295  fprintf(fpl, "%lf %lf %lf %lf %lf %lf\n",
296  beam[0].x(), beam[0].y(), beam[0].z(),
297  beam[nB - 1].x(), beam[nB - 1].y(), beam[nB - 1].z());
298  fprintf(fpv, "%lf %lf %lf\n", beam[0].x(), beam[0].y(), beam[0].z());
299  fprintf(fpv, "%lf %lf %lf\n", beam[nB - 1].x(), beam[nB - 1].y(), beam[nB - 1].z());
300  }
301  }
302  }
303 
304  fclose(fpl);
305  fclose(fpv);
306  fclose(fpm);
307 }
308 
309 /*
310 * GnuPltCubicBentBeam - computes cubic deflection functions from end deflections
311 * and end rotations. Saves deflected shapes to a file. These bent shapes
312 * are exact for mode-shapes, and for frames loaded at their nodes.
313 * Nov/25/2015
314 */
315 void StiffnessIO::GnuPltCubicBentBeam(
316  vector<point> &beam,
317  VX &D,
318  int dual_i,
319  DualGraph *ptr_dualgraph, WireFrame *ptr_frame,
320  double exagg
321  )
322 {
323  double t0, t1, t2, t3, t4, t5, t6, t7, t8, /* coord transf matrix entries */
324  u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12,
325  s, v, w, dX, dY, dZ;
326 
327  WF_edge *ei = ptr_frame->GetEdge(ptr_dualgraph->e_orig_id(dual_i));
328  WF_edge *ej = ei->ppair_;
329  int dual_u = ptr_dualgraph->v_dual_id(ej->pvert_->ID());
330  int dual_v = ptr_dualgraph->v_dual_id(ei->pvert_->ID());
331 
332  double L = ei->Length();
333 
334  int i1, i2;
335  int info;
336  char errMsg[MAXL];
337 
338  MX A(4,4);
339  VX a(4);
340  VX b(4);
341 
342  trsf_.CreateTransMatrix(ej->pvert_->Position(), ei->pvert_->Position(),
343  t0, t1, t2, t3, t4, t5, t6, t7, t8, 0.0);
344 
345  i1 = 6 * dual_u; i2 = 6 * dual_v;
346 
347  /* compute end deflections in local coordinates */
348 
349  u1 = exagg*(t0*D[i1 + 0] + t1*D[i1 + 1] + t2*D[i1 + 2]);
350  u2 = exagg*(t3*D[i1 + 0] + t4*D[i1 + 1] + t5*D[i1 + 2]);
351  u3 = exagg*(t6*D[i1 + 0] + t7*D[i1 + 1] + t8*D[i1 + 2]);
352 
353  u4 = exagg*(t0*D[i1 + 3] + t1*D[i1 + 4] + t2*D[i1 + 5]);
354  u5 = exagg*(t3*D[i1 + 3] + t4*D[i1 + 4] + t5*D[i1 + 5]);
355  u6 = exagg*(t6*D[i1 + 3] + t7*D[i1 + 4] + t8*D[i1 + 5]);
356 
357  u7 = exagg*(t0*D[i2 + 0] + t1*D[i2 + 1] + t2*D[i2 + 2]);
358  u8 = exagg*(t3*D[i2 + 0] + t4*D[i2 + 1] + t5*D[i2 + 2]);
359  u9 = exagg*(t6*D[i2 + 0] + t7*D[i2 + 1] + t8*D[i2 + 2]);
360 
361  u10 = exagg*(t0*D[i2 + 3] + t1*D[i2 + 4] + t2*D[i2 + 5]);
362  u11 = exagg*(t3*D[i2 + 3] + t4*D[i2 + 4] + t5*D[i2 + 5]);
363  u12 = exagg*(t6*D[i2 + 3] + t7*D[i2 + 4] + t8*D[i2 + 5]);
364 
365  /* curve-fitting problem for a cubic polynomial */
366 
367  a[0] = u2; b[0] = u3;
368  a[1] = u8; b[1] = u9;
369  a[2] = u6; b[2] = -u5;
370  a[3] = u12; b[3] = -u11;
371 
372  u7 += L;
373  A(0,0) = 1.0; A(0,1) = u1; A(0,2) = u1*u1; A(0,3) = u1*u1*u1;
374  A(1,0) = 1.0; A(1,1) = u7; A(1,2) = u7*u7; A(1,3) = u7*u7*u7;
375  A(2,0)= 0.0; A(2,1) = 1.; A(2,2) = 2.*u1; A(2,3) = 3.*u1*u1;
376  A(3,0) = 0.0; A(3,1) = 1.; A(3,2) = 2.*u7; A(3,3) = 3.*u7*u7;
377  u7 -= L;
378 
379  VX xa(4);
380  info = solver_.LUDecomp(A, xa, a); /* solve for cubic coef's */
381 
382  if (!info)
383  {
384  sprintf(errMsg, " n1 = %d n2 = %d L = %e u7 = %e \n", dual_u+1, dual_v+1, L, u7);
385  printf(errMsg);
386  exit(30);
387  }
388 
389  VX xb(4);
390  info = solver_.LUDecomp(A, xb, b); /* solve for cubic coef's */
391 
392  if (!info)
393  {
394  sprintf(errMsg, " n1 = %d n2 = %d L = %e u7 = %e \n", dual_u + 1, dual_v + 1, L, u7);
395  printf(errMsg);
396  exit(30);
397  }
398 
399  // debug ... if deformed mesh exageration is too big, some elements
400  // may not be plotted.
401  //fprintf( fpm, "# u1=%e L+u7=%e, dx = %e \n",
402  // u1, fabs(L+u7), fabs(L+u7-u1)/10.0);
403  for (s = u1; fabs(s) <= 1.01*fabs(L + u7); s += fabs(L + u7 - u1) / 10.0)
404  {
405 
406  /* deformed shape in local coordinates */
407  v = xa[0] + xa[1] * s + xa[2] * s*s + xa[3] * s*s*s;
408  w = xb[0] + xb[1] * s + xb[2] * s*s + xb[3] * s*s*s;
409 
410  /* deformed shape in global coordinates */
411  dX = t0*s + t3*v + t6*w;
412  dY = t1*s + t4*v + t7*w;
413  dZ = t2*s + t5*v + t8*w;
414 
415  beam.push_back(point(ej->pvert_->Position().x() + dX,
416  ej->pvert_->Position().y() + dY, ej->pvert_->Position().z() + dZ));
417  //fprintf(fpm, " %12.4e %12.4e %12.4e\n",);
418  }
419 
420  return;
421 }
422 
423 
424 void StiffnessIO::WriteInputData(
425  const char *fpath,
426  DualGraph *ptr_dualgraph,
427  FiberPrintPARM *ptr_parm,
428  int verbose
429  )
430 {
431  FILE *fp;
432  string title_s = "FiberPrint Test File -- static analysis (N,mm,Ton)\n";
433  char errMsg[512];
434 
435  if ((fp = fopen(fpath, "w")) == NULL)
436  {
437  sprintf(errMsg, "\n ERROR: cannot open .3dd transfer data file '%s'", fpath);
438  printf(errMsg);
439  exit(11);
440  }
441 
442  fprintf(fp, title_s.c_str());
443  fprintf(fp, "\n");
444 
445  int nN = ptr_dualgraph->SizeOfFaceList();
446  int nE = ptr_dualgraph->SizeOfVertList();
447  WireFrame *ptr_wf = ptr_dualgraph->ptr_frame_;
448  std::vector<WF_edge*> wf_edge_list = *ptr_wf->GetEdgeList();
449  std::vector<DualFace*> dual_face_list = *ptr_dualgraph->GetFaceList();
450 
451  double r = ptr_parm->radius_;
452  double nr = 0.0;
453  double density = ptr_parm->density_;
454  double g = ptr_parm->g_;
455  double G = ptr_parm->shear_modulus_;
456  double E = ptr_parm->youngs_modulus_;
457  double v = ptr_parm->poisson_ratio_;
458 
459  // Write node data
460  fprintf(fp, "%d # number of nodes\n", nN);
461  fprintf(fp, "#.node x y z r\n");
462  fprintf(fp, "# mm mm mm m\n\n");
463  for (int i = 0; i < nN; i++)
464  {
465  int o_id = ptr_dualgraph->v_orig_id(i);
466 
467  fprintf(fp, "%d %f %f %f %f\n", i + 1,
468  ptr_wf->GetPosition(o_id).x(), ptr_wf->GetPosition(o_id).y(), ptr_wf->GetPosition(o_id).z(),0.0);
469  }
470 
471  // Write nodes with reactions
472  int nR = 0; // Number of restrained nodes
473  std::vector<int> res_index;
474  for (int i = 0; i < nN; i++)
475  {
476  int id = dual_face_list[i]->orig_id();
477 
478  if (ptr_wf->isFixed(id))
479  {
480  nR++;
481  res_index.push_back(i+1);
482  }
483  }
484  fprintf(fp, "\n");
485  fprintf(fp, "%d # number of nodes with reaction\n", nR);
486  fprintf(fp, "#.node x y z xx yy zz 1=b_fixed, 0=free\n");
487  for (int i = 0; i < nR; i++)
488  {
489  fprintf(fp, "%d 1 1 1 1 1 1\n", res_index[i]);
490  }
491 
492  // Write Frame Element data
493  double Ax = F_PI * r * r;
494  double Asy = Ax * (6 + 12 * v + 6 * v*v) / (7 + 12 * v + 4 * v*v);
495  double Asz = Asy;
496  double Jxx = 0.5 * M_PI * r * r * r * r;
497  double Ksy = 0; // shear deformation constant
498  double Ksz = 0;
499 
500  fprintf(fp, "\n");
501  fprintf(fp, "%d # number of frame element\n", nE);
502  fprintf(fp, "#.e n1 n2 Ax Asy Asz Jxx Iyy Izz E G roll density\n");
503  fprintf(fp, "# . . mm^2 mm^2 mm^2 mm^4 mm^4 mm^4 MPa MPa deg T/mm^3\n");
504  for (int i = 0; i < nE; i++)
505  {
506  int e_id = ptr_dualgraph->e_orig_id(i);
507  WF_edge *ei = wf_edge_list[e_id];
508  int u = ei->ppair_->pvert_->ID();
509  int v = ei->pvert_->ID();
510  double L = ei->Length();
511 
512  int dual_u = ptr_dualgraph->v_dual_id(u);
513  int dual_v = ptr_dualgraph->v_dual_id(v);
514 
515  double Iyy = F_PI * r * r * r * r / 4;
516  double Izz = Iyy;
517 
518  fprintf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f %.14f\n",
519  i + 1, dual_u + 1, dual_v + 1,
520  Ax, Asy, Asz, Jxx, Iyy, Izz, E, G, 0.0, density);
521  }
522 
523  printf("\n\n");
524 
525  // parse option for stiffness matrix
526  fprintf(fp, "%d # 1: include shear deformation\n", 0);
527  fprintf(fp, "%d # 1: include geometric stiffness\n", 0);
528  fprintf(fp, "%.1f # exaggerate static mesh deformation\n", 10.0);
529  fprintf(fp, "%.1f # zoom scale for 3D plotting\n", 2.5);
530  fprintf(fp, "%d # x-axis increment for internal forces, mm\n", -1);
531  fprintf(fp, " # if dx is -1 then internal force calculation are skipped\n");
532 
533  // load case parsing
534  fprintf(fp, "\n");
535  fprintf(fp, "%d # number of static load cases\n", 1);
536  fprintf(fp, " # Begin static Load Case 1 of 1\n");
537 
538  fprintf(fp, "\n");
539  fprintf(fp, "# gravitational acceleration for self-weight loading (global)\n");
540  fprintf(fp, "#.gX gY gZ\n");
541  fprintf(fp, "#.mm/s^2 mm/s^2 mm/s^2\n");
542  fprintf(fp, " 0 0 %.6f\n", g);
543 
544  fprintf(fp,"\n");
545  fprintf(fp, "%d # number of loaded nodes\n", 0);
546  fprintf(fp, "%d # number of uniform loads\n", 0);
547  fprintf(fp, "%d # number of trapezoidal loads\n", 0);
548  fprintf(fp, "%d # number of internal concentrated loads\n", 0);
549  fprintf(fp, "%d # number of temperature loads\n", 0);
550 
551  fprintf(fp, "%d # number of nodes with prescribed displacement\n", nR);
552  fprintf(fp, "#.node X-disp1 Y-disp1 Z-disp1 X-rot'n Y-rot'n Z-rot'n\n");
553  fprintf(fp, "# mm mm mm radian radian radian\n");
554  for (int i = 0; i < nR; i++)
555  {
556  fprintf(fp, "%d 0 0 0 0 0 0\n", res_index[i]);
557  }
558 
559  fprintf(fp, " # End static Load Case 1 of 1\n");
560 
561  // Dynamic Modes
562  fprintf(fp, "%d #number of dynamic modes\n",0);
563  fprintf(fp, "# End of transfer data file for fiber test");
564 
565  fclose(fp);
566 }
567 
568 /*
569 * SaveUpperMatrix - save a symmetric matrix of dimension [1..n][1..n] Nov/26/2015
570 * to the named file, use only upper-triangular part
571 */
572 void StiffnessIO::SaveUpperMatrix(char filename[], const MX &A, int n)
573 {
574  FILE *fp_m;
575  int i, j;
576  time_t now;
577 
578  if ((fp_m = fopen(filename, "w")) == NULL)
579  {
580  printf(" error: cannot open file: %s \n", filename);
581  exit(1016);
582  }
583 
584  (void)time(&now);
585  fprintf(fp_m, "%% filename: %s - %s\n", filename, ctime(&now));
586  fprintf(fp_m, "%% type: matrix \n");
587  fprintf(fp_m, "%% rows: %d\n", n);
588  fprintf(fp_m, "%% columns: %d\n", n);
589  for (i = 0; i < n; i++)
590  {
591  for (j = 0; j < n; j++)
592  {
593  if (i > j)
594  {
595  if (fabs(A(j,i)) > 1.e-99)
596  {
597  fprintf(fp_m, "%21.12e", A(j,i));
598  }
599  else
600  {
601  fprintf(fp_m, " 0 ");
602  }
603  }
604  else
605  {
606  if (fabs(A(i, j)) > 1.e-99)
607  {
608  fprintf(fp_m, "%21.12e", A(i,j));
609  }
610  else
611  {
612  fprintf(fp_m, " 0 ");
613  }
614  }
615  }
616  fprintf(fp_m, "\n");
617  }
618  fclose(fp_m);
619  return;
620 
621 }
622 
623 void StiffnessIO::SaveDisplaceVector(char filename[], const VX &D, int n, DualGraph *ptr_dual_graph)
624 {
625  FILE *fp;
626  int i, j;
627  time_t now;
628  int nN = n / 6;
629  vector<DualFace*> dual_face_list = *ptr_dual_graph->GetFaceList();
630  VX tmp_D(D.size());
631 
632  if ((fp = fopen(filename, "w")) == NULL)
633  {
634  printf(" error: cannot open file: %s \n", filename);
635  exit(1016);
636  }
637 
638  (void)time(&now);
639  fprintf(fp, "%% filename: %s - %s\n", filename, ctime(&now));
640  fprintf(fp, "%% type: vector \n");
641  fprintf(fp, "%% rows: %d\n", n);
642 
643  fprintf(fp, "\n");
644  fprintf(fp, "NODE DISPLACEMENTS (global)\n");
645  fprintf(fp, "#.node X-dsp Y-dsp Z-dsp X-rot Y-rot Z-rot\n");
646  for (i = 0; i < nN; i++)
647  {
648  /* i is the dual face id, convert it back to wf_vert id*/
649  int v_id = dual_face_list[i]->orig_id();
650  int v_dual_id = ptr_dual_graph->v_dual_id(v_id);
651  for (j = 0; j < 6; j++)
652  {
653  tmp_D[6 * v_dual_id + j] = D[6 * i + j];
654  }
655  }
656 
657  for (i = 0; i < nN; i++)
658  {
659  fprintf(fp, "%d %.6f %.6f %.6f %.6f %.6f %.6f\n",
660  i+1, tmp_D[6 * i], tmp_D[6 * i + 1], tmp_D[6 * i + 2], tmp_D[6 * i + 3], tmp_D[6 * i + 4], tmp_D[6 * i + 5]);
661  }
662 
663  fclose(fp);
664  return;
665 
666 }
u2
GLfixed GLfixed u2
n
GLdouble n
DualGraph::GetFaceList
vector< DualFace * > * GetFaceList()
Definition: DualGraph.h:143
DualGraph::v_orig_id
int v_orig_id(int i)
Definition: DualGraph.h:154
StiffnessIO::GnuPltCubicBentBeam
void GnuPltCubicBentBeam(vector< point > &beam, VX &D, int dual_i, DualGraph *ptr_dualgraph, WireFrame *ptr_frame, double exagg)
Definition: StiffnessIO.cpp:315
u1
GLfixed u1
m
const GLfloat * m
point
trimesh::point point
Definition: WireFrame.h:51
t0
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
scale
GLenum GLenum GLenum GLenum GLenum scale
void
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
StiffnessIO::trsf_
CoordTrans trsf_
Definition: StiffnessIO.h:159
t1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
WireFrame::GetEdge
WF_edge * GetEdge(int i)
Definition: WireFrame.h:238
StiffnessIO::GnuPltStaticMesh
void GnuPltStaticMesh(const char *fpath, const char *meshpath, const char *plotpath, VX &D, double exagg_static, float scale, DualGraph *ptr_dualgraph, WireFrame *ptr_frame)
Definition: StiffnessIO.cpp:46
trimesh::Vec::x
reference x()
Definition: Vec.h:194
StiffnessIO::VX
Eigen::VectorXd VX
Definition: StiffnessIO.h:74
WF_edge
Definition: WireFrame.h:113
DualGraph::SizeOfFaceList
int SizeOfFaceList()
Definition: DualGraph.h:147
M_PI
#define M_PI
y2
GLfixed GLfixed GLfixed y2
WireFrame::GetEdgeList
vector< WF_edge * > * GetEdgeList()
Definition: WireFrame.h:236
x
GLint GLenum GLint x
len
GLenum GLsizei len
w
GLubyte GLubyte GLubyte GLubyte w
a
GLboolean GLboolean GLboolean GLboolean a
MAXL
#define MAXL
Definition: GCommon.h:44
WF_edge::Length
double Length() const
Definition: WireFrame.h:142
g
GLboolean GLboolean g
DualGraph::e_orig_id
int e_orig_id(int u)
Definition: DualGraph.h:152
WF_edge::ppair_
WF_edge * ppair_
Definition: WireFrame.h:166
y1
GLfixed y1
VERSION
#define VERSION
Definition: GCommon.h:50
FILENMAX
#define FILENMAX
Definition: GCommon.h:47
StiffnessIO::WriteInputData
void WriteInputData(const char *fpath, DualGraph *ptr_dualgraph, FiberPrintPARM *ptr_parm, int verbose)
Definition: StiffnessIO.cpp:424
StiffnessIO::MX
Eigen::MatrixXd MX
Definition: StiffnessIO.h:72
F_PI
#define F_PI
Definition: GCommon.h:54
FiberPrintPARM::youngs_modulus_
double youngs_modulus_
Definition: FiberPrintPARM.h:71
nx
GLbyte nx
FiberPrintPARM::shear_modulus_
double shear_modulus_
Definition: FiberPrintPARM.h:72
trimesh::Vec::y
reference y()
Definition: Vec.h:203
FiberPrintPARM::poisson_ratio_
double poisson_ratio_
Definition: FiberPrintPARM.h:73
StiffnessIO::SaveUpperMatrix
void SaveUpperMatrix(char filename[], const MX &A, int n)
Definition: StiffnessIO.cpp:572
b
GLboolean GLboolean GLboolean b
z
GLdouble GLdouble GLdouble z
DualGraph::ptr_frame_
WireFrame * ptr_frame_
Definition: DualGraph.h:189
WF_vert::ID
int ID() const
Definition: WireFrame.h:80
x1
GLuint GLfloat GLfloat GLfloat x1
CoordTrans::CreateTransMatrix
void CreateTransMatrix(std::vector< V3 > xyz, double L, int n1, int n2, double &t0, double &t1, double &t2, double &t3, double &t4, double &t5, double &t6, double &t7, double &t8, float p)
Definition: CoordTrans.cpp:3
v
const GLdouble * v
r
GLboolean r
WF_edge::pvert_
WF_vert * pvert_
Definition: WireFrame.h:164
s
GLdouble s
WireFrame::isFixed
bool isFixed(int u) const
Definition: WireFrame.h:253
StiffnessIO::solver_
StiffnessSolver solver_
Definition: StiffnessIO.h:160
x2
GLfixed GLfixed x2
DualGraph::v_dual_id
int v_dual_id(int i)
Definition: DualGraph.h:155
StiffnessIO::SaveDisplaceVector
void SaveDisplaceVector(char filename[], const VX &D, int n, DualGraph *ptr_dual_graph)
Definition: StiffnessIO.cpp:623
DualGraph::SizeOfVertList
int SizeOfVertList()
Definition: DualGraph.h:145
StiffnessIO::OutputPath
void OutputPath(const char *fname, char fullpath[], const int len, char *default_outdir, int verbose)
Definition: StiffnessIO.cpp:7
res
GLuint res
trimesh::Vec::z
reference z()
Definition: Vec.h:212
WireFrame::GetPosition
point GetPosition(int u) const
Definition: WireFrame.h:241
y
GLint y
WF_vert::Position
point Position() const
Definition: WireFrame.h:78
StiffnessSolver::LUDecomp
bool LUDecomp(MX &A, VX &x, VX &b)
Definition: StiffnessSolver.cpp:135
StiffnessIO::TempDir
const char * TempDir()
Definition: StiffnessIO.h:144

choreo_task_sequence_planner
Author(s): Yijiang Huang
autogenerated on Thu Jul 18 2019 04:03:14