rsvs3D  0.0.0
Codes for the c++ implementation of the 3D RSVS
postprocessing.cpp
1 #include "postprocessing.hpp"
2 
3 #include "RSVScalc.hpp"
4 #include "mesh.hpp"
5 #include "meshprocessing.hpp"
6 #include "triangulate.hpp"
7 #include "warning.hpp"
8 #include <iostream>
9 #include <vector>
10 
11 namespace tecplotconst = rsvs3d::constants::tecplot;
12 
13 using namespace std;
14 
15 // Functions
16 void ExtractMeshData(const mesh &grid, int *nVert, int *nEdge, int *nVolu, int *nSurf, int *totNumFaceNode)
17 {
18  // Extracts Data needed to write out a mesh to tecplot
19  int ii;
20 
21  *nVert = grid.verts.size();
22  *nVolu = grid.volus.size();
23  *nSurf = grid.surfs.size();
24  *nEdge = grid.edges.size();
25  *totNumFaceNode = 0;
26  for (ii = 0; ii < *nSurf; ++ii)
27  {
28  *totNumFaceNode = *totNumFaceNode + int(grid.surfs(ii)->edgeind.size());
29  }
30 }
31 
32 namespace dataoutput
33 {
34 void Coord(tecplotfile &tecout, const mesh &meshout, int nVert, int nVertDat)
35 {
36  int ii, jj, nCoord;
37  // Print vertex Data
38  nCoord = int(meshout.verts(0)->coord.size());
39  nCoord = nCoord > nVertDat ? nVertDat : nCoord;
40  for (jj = 0; jj < nCoord; ++jj)
41  {
42  for (ii = 0; ii < nVert; ++ii)
43  {
44  tecout.Print("%.16lf ", meshout.verts(ii)->coord[jj]);
45  }
46  tecout.NewLine();
47  }
48  for (jj = int(meshout.verts(0)->coord.size()); jj < nVertDat; ++jj)
49  {
50  for (ii = 0; ii < nVert; ++ii)
51  {
52  tecout.Print("%lf ", 0.0);
53  }
54  tecout.NewLine();
55  }
56 }
57 
58 void Snaxel(tecplotfile &tecout, const snake &snakeout, int nVert)
59 {
60  int ii;
61  for (ii = 0; ii < nVert; ++ii)
62  {
63  tecout.Print("%.16lf ", snakeout.snaxs(ii)->d);
64  }
65  tecout.NewLine();
66  for (ii = 0; ii < nVert; ++ii)
67  {
68  tecout.Print("%.16lf ", snakeout.snaxs(ii)->v);
69  }
70  tecout.NewLine();
71  for (ii = 0; ii < nVert; ++ii)
72  {
73  auto temp = snakeout.snakeconn.verts(ii)->elmind(snakeout.snakeconn, 2);
74  int maxSize = 0;
75  for (auto surfSub : temp)
76  {
77  int tempSize = snakeout.snakeconn.surfs.isearch(surfSub)->edgeind.size();
78  maxSize = maxSize < tempSize ? tempSize : maxSize;
79  }
80  // tecout.Print("%i ",snakeout.snaxs(ii)->isfreeze);
81  tecout.Print("%i ", maxSize);
82  }
83  tecout.NewLine();
84 }
85 
93 void VertexNormal(tecplotfile &tecout, const mesh &meshin, int nVert)
94 {
95  std::vector<double> coords = MeshUnitNormals(meshin);
96 
97  for (int j = 0; j < 3; ++j)
98  {
99  for (int i = 0; i < nVert; ++i)
100  {
101  tecout.Print("%.16lf ", coords[i * 3 + j]);
102  }
103  tecout.NewLine();
104  }
105 }
106 void VertexLaplacian(tecplotfile &tecout, const mesh &meshin, int nVert)
107 {
108  auto coords = MeshLaplacians(meshin);
109  for (int j = 0; j < 3; ++j)
110  {
111  for (int i = 0; i < nVert; ++i)
112  {
113  tecout.Print("%.16lf ", coords[i * 3 + j]);
114  }
115  tecout.NewLine();
116  }
117 }
118 
119 void SnaxelDirection(tecplotfile &tecout, const snake &snakein, int nVert)
120 {
121  auto coords = snakein.MoveDirections();
122  for (int j = 0; j < 3; ++j)
123  {
124  for (int i = 0; i < nVert; ++i)
125  {
126  tecout.Print("%.16lf ", coords[i * 3 + j]);
127  }
128  tecout.NewLine();
129  }
130 }
131 
132 } // namespace dataoutput
133 
134 int tecplotfile::VolDataBlock(const mesh &meshout, int nVert, int nVolu, int nVertDat, const std::vector<int> &voluList,
135  const std::vector<int> &vertList)
136 {
137  // Prints the Coord and Fill Data blocks to the tecplot file
138 
139  int ii, jj, nCoord;
140  // Print vertex Data
141  nCoord = int(meshout.verts(0)->coord.size());
142  nCoord = nCoord > nVertDat ? nVertDat : nCoord;
143  if (vertList.size() == 0)
144  {
145  dataoutput::Coord(*this, meshout, nVert, nVertDat);
146  }
147  else
148  {
149  for (jj = 0; jj < nCoord; ++jj)
150  {
151  for (auto ii : vertList)
152  {
153  this->Print("%.16lf ", meshout.verts.isearch(ii)->coord[jj]);
154  }
155  this->NewLine();
156  }
157  nVert = vertList.size();
158  for (jj = int(meshout.verts(0)->coord.size()); jj < nVertDat; ++jj)
159  {
160  for (ii = 0; ii < nVert; ++ii)
161  {
162  this->Print("%lf ", 0.0);
163  }
164  this->NewLine();
165  }
166  }
167  // Print Cell Data
168  if (voluList.size() == 0)
169  {
170  for (ii = 0; ii < nVolu; ++ii)
171  {
172  this->Print("%.16lf ", meshout.volus(ii)->fill);
173  }
174  this->NewLine();
175  for (ii = 0; ii < nVolu; ++ii)
176  {
177  this->Print("%.16lf ", meshout.volus(ii)->target);
178  }
179  this->NewLine();
180  for (ii = 0; ii < nVolu; ++ii)
181  {
182  this->Print("%.16lf ", meshout.volus(ii)->error);
183  }
184  this->NewLine();
185  }
186  else
187  {
188  for (auto ii : voluList)
189  {
190  this->Print("%.16lf ", meshout.volus.isearch(ii)->fill);
191  }
192  this->NewLine();
193 
194  for (auto ii : voluList)
195  {
196  this->Print("%.16lf ", meshout.volus.isearch(ii)->target);
197  }
198  this->NewLine();
199  for (auto ii : voluList)
200  {
201  this->Print("%.16lf ", meshout.volus.isearch(ii)->error);
202  }
203  this->NewLine();
204  }
205  return (0);
206 }
207 
208 int tecplotfile::SnakeDataBlock(const snake &snakeout, int nVert, int nVertDat, std::string snakeData, bool printCoord)
209 {
210  // Prints the Coord and Fill Data blocks to the tecplot file
211  if (printCoord)
212  {
213  dataoutput::Coord(*this, snakeout.snakeconn, nVert, nVertDat);
214  }
215  // Print Cell Data
216  if (snakeData.compare(tecplotconst::snakedata::snaxel) == 0)
217  {
218  dataoutput::Snaxel(*this, snakeout, nVert);
219  }
220  else if (snakeData.compare(tecplotconst::snakedata::normal) == 0)
221  {
222  dataoutput::VertexNormal(*this, snakeout.snakeconn, nVert);
223  }
224  else if (snakeData.compare(tecplotconst::snakedata::laplacian) == 0)
225  {
226  dataoutput::VertexLaplacian(*this, snakeout.snakeconn, nVert);
227  }
228  else if (snakeData.compare(tecplotconst::snakedata::direction) == 0)
229  {
230  dataoutput::SnaxelDirection(*this, snakeout, nVert);
231  }
232  else
233  {
234  stringstream errstr;
235  errstr << "Unknown snake data output '" << snakeData << "'";
236  RSVS3D_ERROR_ARGUMENT(errstr.str().c_str());
237  }
238 
239  return (0);
240 }
241 
242 int tecplotfile::SurfDataBlock(const mesh &meshout, int nVert, int nSurf, int nVertDat)
243 {
244  // Prints the Coord and Fill Data blocks to the tecplot file
245 
246  int ii;
247  // Print vertex Data
248  dataoutput::Coord(*this, meshout, nVert, nVertDat);
249  // Print Cell Data
250  for (ii = 0; ii < nSurf; ++ii)
251  {
252  this->Print("%.16lf ", meshout.surfs(ii)->fill);
253  }
254  this->NewLine();
255  for (ii = 0; ii < nSurf; ++ii)
256  {
257  this->Print("%.16lf ", meshout.surfs(ii)->target);
258  }
259  this->NewLine();
260  for (ii = 0; ii < nSurf; ++ii)
261  {
262  this->Print("%.16lf ", meshout.surfs(ii)->error);
263  }
264  this->NewLine();
265  return (0);
266 }
267 
268 int tecplotfile::LineDataBlock(const mesh &meshout, int nVert, int nEdge, int nVertDat, int nCellDat)
269 {
270  // Prints the Coord and Fill Data blocks to the tecplot file
271 
272  int ii, jj;
273  // Print vertex Data
274  dataoutput::Coord(*this, meshout, nVert, nVertDat);
275  // Print Cell Data
276  for (jj = 0; jj < nCellDat; ++jj)
277  {
278  for (ii = 0; ii < nEdge; ++ii)
279  {
280  this->Print("%i ", 0);
281  }
282  this->NewLine();
283  }
284  return (0);
285 }
286 
287 int tecplotfile::VertDataBlock(const mesh &meshout, int nVert, int nVertDat, int nCellDat, const vector<int> &vertList)
288 {
289  // Prints the Coord and Fill Data blocks to the tecplot file
290 
291  int ii, jj, nCoord;
292  nCoord = int(meshout.verts(0)->coord.size());
293  // Print vertex Data
294  if (vertList.size() > 0)
295  { // if vertList is a list of index
296  for (jj = 0; jj < nCoord; ++jj)
297  {
298  for (ii = 0; ii < nVert; ++ii)
299  {
300  this->Print("%.16lf ", meshout.verts.isearch(vertList[ii])->coord[jj]);
301  }
302  this->NewLine();
303  }
304  }
305  else if (int(vertList.size()) == meshout.verts.size())
306  { // vertList is a boolean
307  for (jj = 0; jj < nCoord; ++jj)
308  {
309  for (ii = 0; ii < nVert; ++ii)
310  {
311  if (vertList[ii])
312  {
313  this->Print("%.16lf ", meshout.verts(ii)->coord[jj]);
314  }
315  }
316  this->NewLine();
317  }
318  }
319  else
320  {
321  for (jj = 0; jj < nCoord; ++jj)
322  {
323  for (ii = 0; ii < nVert; ++ii)
324  {
325  this->Print("%.16lf ", meshout.verts(ii)->coord[jj]);
326  }
327  this->NewLine();
328  }
329  }
330  for (jj = nCoord; jj < nVertDat; ++jj)
331  {
332  for (ii = 0; ii < nVert; ++ii)
333  {
334  this->Print("%lf ", 0.0);
335  }
336  this->NewLine();
337  }
338  // Print Cell Data
339  for (jj = 0; jj < nCellDat; ++jj)
340  {
341  for (ii = 0; ii < nVert; ++ii)
342  {
343  this->Print("%i ", 0);
344  }
345  this->NewLine();
346  }
347  return (0);
348 }
349 
350 int tecplotfile::SurfFaceMap(const mesh &meshout, int nEdge)
351 {
352  int ii, jj, actVert;
353  int verts[2];
354 
355  for (ii = 0; ii < nEdge; ++ii)
356  { // Print Number of vertices per face
357  verts[0] = meshout.edges(ii)->vertind[0];
358  verts[1] = meshout.edges(ii)->vertind[1];
359  this->Print("%i %i\n", meshout.verts.find(verts[0]) + 1, meshout.verts.find(verts[1]) + 1);
360  this->ResetLine();
361  }
362 
363  for (jj = 0; jj < 2; ++jj)
364  { // print index of left and right facing volumes
365  for (ii = 0; ii < nEdge; ++ii)
366  {
367  // cout << ii << "," << jj << endl ;
368  actVert = meshout.edges(ii)->surfind[jj];
369  this->Print("%i ", meshout.surfs.find(actVert) + 1);
370  }
371  this->NewLine();
372  }
373 
374  return (0);
375 }
376 
377 int tecplotfile::LineFaceMap(const mesh &meshout, int nEdge)
378 {
379  int ii;
380  int verts[2];
381 
382  for (ii = 0; ii < nEdge; ++ii)
383  { // Print Number of vertices per face
384  verts[0] = meshout.edges(ii)->vertind[0];
385  verts[1] = meshout.edges(ii)->vertind[1];
386  this->Print("%i %i\n", meshout.verts.find(verts[0]) + 1, meshout.verts.find(verts[1]) + 1);
387  this->ResetLine();
388  }
389 
390  return (0);
391 }
392 
393 int tecplotfile::VolFaceMap(const mesh &meshout, int nSurf)
394 {
395  int ii, jj, actVert, edgeCurr;
396  int verts[2], vertsPast[2];
397  for (ii = 0; ii < nSurf; ++ii)
398  { // Print Number of vertices per face
399  this->Print("%i ", int(meshout.surfs(ii)->edgeind.size()));
400  }
401  this->NewLine();
402 
403  for (ii = 0; ii < nSurf; ++ii)
404  { // print ordered list of vertices in face
405  jj = int(meshout.surfs(ii)->edgeind.size()) - 1;
406 
407  edgeCurr = meshout.edges.find(meshout.surfs(ii)->edgeind[jj]);
408  verts[0] = meshout.verts.find(meshout.edges(edgeCurr)->vertind[0]);
409  verts[1] = meshout.verts.find(meshout.edges(edgeCurr)->vertind[1]);
410  vertsPast[0] = verts[0];
411  vertsPast[1] = verts[1];
412 
413  for (jj = 0; jj < int(meshout.surfs(ii)->edgeind.size()); ++jj)
414  {
415  edgeCurr = meshout.edges.find(meshout.surfs(ii)->edgeind[jj]);
416 
417  verts[0] = meshout.verts.find(meshout.edges(edgeCurr)->vertind[0]);
418  verts[1] = meshout.verts.find(meshout.edges(edgeCurr)->vertind[1]);
419 
420  if ((verts[0] == vertsPast[0]) || (verts[1] == vertsPast[0]))
421  {
422  actVert = 0;
423  }
424 #ifdef TEST_POSTPROCESSING
425  else if ((verts[0] == vertsPast[1]) || (verts[1] == vertsPast[1]))
426  {
427  actVert = 1;
428  }
429 #endif // TEST_POSTPROCESSING
430  else
431  {
432  actVert = 1;
433 #ifdef TEST_POSTPROCESSING
434  // meshout.surfs(ii)->disptree(meshout,4);
435 
436  cerr << "Warning: postprocessing.cpp:tecplotfile::VolFaceMap" << endl;
437  cerr << " Mesh Output failed in output of facemap data" << endl;
438  cerr << " Surface is not ordered " << endl;
439  return (-1);
440 #endif
441  }
442 
443  this->Print("%i ", vertsPast[actVert] + 1);
444  vertsPast[0] = verts[0];
445  vertsPast[1] = verts[1];
446  }
447  this->NewLine();
448  }
449  for (jj = 0; jj < 2; ++jj)
450  { // print index of left and right facing volumes
451  for (ii = 0; ii < nSurf; ++ii)
452  {
453  // cout << ii << "," << jj << endl ;
454  actVert = meshout.surfs(ii)->voluind[jj];
455  this->Print("%i ", meshout.volus.find(actVert) + 1);
456  }
457  this->NewLine();
458  }
459 
460  return (0);
461 }
462 int tecplotfile::VolFaceMap(const mesh &meshout, const std::vector<int> &surfList, const std::vector<int> &voluList,
463  const std::vector<int> &vertList)
464 {
465  int jj, actVert, edgeCurr;
466  int verts[2], vertsPast[2];
467  for (auto ii : meshout.surfs.find_list(surfList))
468  { // Print Number of vertices per face
469  this->Print("%i ", int(meshout.surfs(ii)->edgeind.size()));
470  }
471  this->NewLine();
472  HashedVector<int, int> verthash;
473  verthash.vec = vertList;
474  verthash.GenerateHash();
475  for (auto ii : meshout.surfs.find_list(surfList))
476  { // print ordered list of vertices in face
477  jj = int(meshout.surfs(ii)->edgeind.size()) - 1;
478 
479  edgeCurr = meshout.edges.find(meshout.surfs(ii)->edgeind[jj]);
480  verts[0] = verthash.find(meshout.edges(edgeCurr)->vertind[0]);
481  verts[1] = verthash.find(meshout.edges(edgeCurr)->vertind[1]);
482  vertsPast[0] = verts[0];
483  vertsPast[1] = verts[1];
484 
485  for (jj = 0; jj < int(meshout.surfs(ii)->edgeind.size()); ++jj)
486  {
487  edgeCurr = meshout.edges.find(meshout.surfs(ii)->edgeind[jj]);
488 
489  verts[0] = verthash.find(meshout.edges(edgeCurr)->vertind[0]);
490  verts[1] = verthash.find(meshout.edges(edgeCurr)->vertind[1]);
491 
492  if ((verts[0] == vertsPast[0]) || (verts[1] == vertsPast[0]))
493  {
494  actVert = 0;
495  }
496 #ifdef TEST_POSTPROCESSING
497  else if ((verts[0] == vertsPast[1]) || (verts[1] == vertsPast[1]))
498  {
499  actVert = 1;
500  }
501 #endif // TEST_POSTPROCESSING
502  else
503  {
504  actVert = 1;
505 #ifdef TEST_POSTPROCESSING
506  // meshout.surfs(ii)->disptree(meshout,4);
507 
508  cerr << "Warning: postprocessing.cpp:tecplotfile::VolFaceMap" << endl;
509  cerr << " Mesh Output failed in output of facemap data" << endl;
510  cerr << " Surface is not ordered " << endl;
511  return (-1);
512 #endif
513  }
514 
515  this->Print("%i ", vertsPast[actVert] + 1);
516  vertsPast[0] = verts[0];
517  vertsPast[1] = verts[1];
518  }
519  this->NewLine();
520  }
521  for (jj = 0; jj < 2; ++jj)
522  { // print index of left and right facing volumes
523  for (auto ii : meshout.surfs.find_list(surfList))
524  {
525  // cout << ii << "," << jj << endl ;
526  actVert = meshout.surfs(ii)->voluind[jj];
527  int k = 0;
528  bool notFoundCell = true;
529  for (auto j : voluList)
530  {
531  k++;
532  if (j == actVert)
533  {
534  notFoundCell = false;
535  break;
536  }
537  }
538  if (notFoundCell)
539  {
540  k = 0;
541  }
542  this->Print("%i ", k);
543  }
544  this->NewLine();
545  }
546 
547  return (0);
548 }
549 // Class function Implementation
550 int tecplotfile::PrintMesh(const mesh &meshout, int strandID, double timeStep, int forceOutType,
551  const vector<int> &vertList)
552 {
553  if (!this->isOpen())
554  {
555  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
556  return -1;
557  }
558  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat;
559 
560  ExtractMeshData(meshout, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
561  if (nVert == 0)
562  { // Don't print mesh
563  return (1);
564  }
565  if (nZones == 0)
566  {
567  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" "
568  ",\"v2\", \"v3\"\n");
569  }
570 
571  this->NewZone();
572 
573  if (strandID > 0)
574  {
575  this->StrandTime(strandID, timeStep);
576  }
577  // Fixed by the dimensionality of the mesh
578  nVertDat = 3;
579  nCellDat = 3;
580 
581  if (forceOutType == tecplotconst::autoselect)
582  {
583  if (nVolu > 0)
584  {
585  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
586  }
587  else if (nSurf > 0)
588  {
589  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
590  }
591  else if (nEdge > 0)
592  {
593  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
594  }
595  else
596  {
597  forceOutType = tecplotconst::point;
598  }
599  }
600 
601  if (forceOutType == tecplotconst::polyhedron)
602  {
603  this->ZoneHeaderPolyhedron(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat);
604  this->VolDataBlock(meshout, nVert, nVolu, nVertDat);
605  this->VolFaceMap(meshout, nSurf);
606  }
607  else if (forceOutType == tecplotconst::polygon)
608  {
609  this->ZoneHeaderPolygon(nVert, nEdge, nSurf, nVertDat, nCellDat);
610  this->SurfDataBlock(meshout, nVert, nSurf, nVertDat);
611  this->SurfFaceMap(meshout, nEdge);
612  }
613  else if (forceOutType == tecplotconst::line)
614  {
615  this->ZoneHeaderFelineseg(nVert, nEdge, nVertDat, nCellDat);
616  this->LineDataBlock(meshout, nVert, nEdge, nVertDat, nCellDat);
617  this->LineFaceMap(meshout, nEdge);
618  }
619  else if (forceOutType == tecplotconst::point)
620  {
621  if (int(vertList.size()) == nVert)
622  {
623  nVert = 0;
624  for (int ii = 0; ii < int(vertList.size()); ++ii)
625  {
626  nVert += int(vertList[ii]);
627  }
628  }
629  else if (vertList.size() > 0)
630  {
631  nVert = int(vertList.size());
632  }
633  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat);
634  this->VertDataBlock(meshout, nVert, nVertDat, nCellDat, vertList);
635  // No map just points
636  }
637  else if (forceOutType == 5)
638  {
639  if (vertList.size() > 0)
640  {
641  nVolu = int(vertList.size());
642  }
643  else
644  {
645  RSVS3D_ERROR_ARGUMENT("voluList with forceOutType == 5 must at"
646  " least be length 1");
647  }
648  nSurf = 0;
649  totNumFaceNode = 0;
650  auto voluSubList = meshout.volus.find_list(vertList);
651  auto surfList = ConcatenateVectorField(meshout.volus, &volu::surfind, voluSubList);
652  sort(surfList);
653  unique(surfList);
654  nSurf = surfList.size();
655  auto surfSubList = meshout.surfs.find_list(surfList);
656  auto edgeList = ConcatenateVectorField(meshout.surfs, &surf::edgeind, surfSubList);
657  auto edgeSubList = meshout.edges.find_list(edgeList);
658  auto actualVert = ConcatenateVectorField(meshout.edges, &edge::vertind, edgeSubList);
659  sort(actualVert);
660  unique(actualVert);
661  nVert = actualVert.size();
662 
663  for (auto j : surfList)
664  {
665  totNumFaceNode += meshout.surfs.isearch(j)->edgeind.size();
666  }
667 
668  this->ZoneHeaderPolyhedron(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat);
669  this->VolDataBlock(meshout, nVert, nVolu, nVertDat, vertList, actualVert);
670  this->VolFaceMap(meshout, surfList, vertList, actualVert);
671  }
672  return (0);
673 }
674 
675 // Class function Implementation
676 int tecplotfile::PrintSnake(const snake &snakeout, int strandID, double timeStep, int forceOutType,
677  const vector<int> &vertList)
678 {
679  return this->PrintSnake(tecplotconst::snakedata::__default, snakeout, strandID, timeStep, forceOutType,
680  tecplotconst::nosharedzone, vertList);
681 }
682 int tecplotfile::PrintSnake(std::string snakeData, const snake &snakeout, int strandID, double timeStep,
683  int forceOutType, int coordConnShareZone, const vector<int> &vertList)
684 {
685  if (!this->isOpen())
686  {
687  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
688  return -1;
689  }
690  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat;
691 
692  ExtractMeshData(snakeout.snakeconn, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
693  if (nVert == 0)
694  { // Don't print mesh
695  return (1);
696  }
697 
698  if (nZones == 0)
699  {
700  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"\n");
701  }
702 
703  this->NewZone();
704 
705  if (strandID > 0)
706  {
707  this->StrandTime(strandID, timeStep);
708  }
709  // Fixed by the dimensionality of the mesh
710  nVertDat = 3;
711  nCellDat = 3;
712 
713  if (forceOutType == tecplotconst::autoselect)
714  {
715  if (nVolu > 0)
716  {
717  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
718  }
719  else if (nSurf > 0)
720  {
721  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
722  }
723  else if (nEdge > 0)
724  {
725  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
726  }
727  else
728  {
729  forceOutType = tecplotconst::point;
730  }
731  }
732 
733  bool zoneSharing = tecplotconst::__issharedzone(coordConnShareZone);
734  bool printCoord = !zoneSharing;
735  if (forceOutType == tecplotconst::polyhedron)
736  {
737  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat);
738  if (zoneSharing)
739  {
740  this->DefShareZoneVolume(coordConnShareZone, nVertDat);
741  }
742  this->SnakeDataBlock(snakeout, nVert, nVertDat, snakeData, printCoord);
743  if (printCoord)
744  {
745  this->VolFaceMap(snakeout.snakeconn, nSurf);
746  }
747  }
748  else if (forceOutType == tecplotconst::polygon)
749  {
750  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat);
751  if (zoneSharing)
752  {
753  this->DefShareZoneVolume(coordConnShareZone, nVertDat);
754  }
755  this->SnakeDataBlock(snakeout, nVert, nVertDat, snakeData, printCoord);
756  if (printCoord)
757  {
758  this->SurfFaceMap(snakeout.snakeconn, nEdge);
759  }
760  }
761  else if (forceOutType == tecplotconst::line)
762  {
763  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat);
764  if (zoneSharing)
765  {
766  this->DefShareZoneVolume(coordConnShareZone, nVertDat);
767  }
768  this->SnakeDataBlock(snakeout, nVert, nVertDat, snakeData, printCoord);
769  if (printCoord)
770  {
771  this->LineFaceMap(snakeout.snakeconn, nEdge);
772  }
773  }
774  else if (forceOutType == tecplotconst::point)
775  {
776  if (int(vertList.size()) == nVert)
777  {
778  nVert = 0;
779  for (int ii = 0; ii < int(vertList.size()); ++ii)
780  {
781  nVert += int(vertList[ii]);
782  }
783  }
784  else if (vertList.size() > 0)
785  {
786  nVert = int(vertList.size());
787  }
788  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat);
789  if (zoneSharing)
790  {
791  this->DefShareZoneVolume(coordConnShareZone, nVertDat);
792  }
793  this->SnakeDataBlock(snakeout, nVert, nVertDat, snakeData, printCoord);
794  // No map just points
795  }
796  return (0);
797 }
798 
799 int tecplotfile::PrintVolumeDat(const mesh &meshout, int shareZone, int strandID, double timeStep)
800 {
801  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat;
802  int forceOutType;
803  if (nZones == 0)
804  {
805  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"\n");
806  }
807 
808  this->NewZone();
809 
810  if (strandID > 0)
811  {
812  this->StrandTime(strandID, timeStep);
813  }
814  ExtractMeshData(meshout, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
815  // Fixed by the dimensionality of the mesh
816  nVertDat = 3;
817  nCellDat = 3;
818  // forceOutType=tecplotconst::autoselect;
819  // if (forceOutType==tecplotconst::autoselect){
820  if (nVolu > 0)
821  {
822  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
823  }
824  else if (nSurf > 0)
825  {
826  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
827  }
828  else if (nEdge > 0)
829  {
830  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
831  }
832  else
833  {
834  forceOutType = tecplotconst::point;
835  }
836  // }
837 
838  if (forceOutType == tecplotconst::polyhedron)
839  {
840  this->ZoneHeaderPolyhedron(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat);
841  this->DefShareZoneVolume(shareZone, nVertDat);
842  this->VolDataBlock(meshout, nVert, nVolu, 0);
843  // this->VolFaceMap(meshout,nSurf);
844  }
845  else if (forceOutType == tecplotconst::polygon)
846  {
847  this->ZoneHeaderPolygon(nVert, nEdge, nSurf, nVertDat, nCellDat);
848  this->DefShareZoneVolume(shareZone, nVertDat);
849  this->SurfDataBlock(meshout, nVert, nSurf, 0);
850  // this->SurfFaceMap(meshout,nEdge);
851  }
852  else if (forceOutType == tecplotconst::line)
853  {
854  RSVS3D_ERROR_ARGUMENT("Cannot output volume of line");
855  }
856  else if (forceOutType == tecplotconst::point)
857  {
858  RSVS3D_ERROR_ARGUMENT("Cannot output volume of point");
859  }
860  return (0);
861  // CONNECTIVITYSHAREZONE=<zone>
862  // VARSHARELIST=([set-of-vars]=<zone>, [set-ofvars]=<zone>)
863 }
864 int tecplotfile::DefShareZoneVolume(int shareZone, int nVertDat)
865 {
866  this->Print("CONNECTIVITYSHAREZONE=%i\n", shareZone);
867  this->Print("VARSHARELIST=([%i-%i]=%i )\n", 1, nVertDat, shareZone);
868 
869  return (0);
870 }
871 
872 int tecplotfile::PrintSnakeInternalPts(const snake &snakein, int strandID, double timeStep)
873 {
874  if (!this->isOpen())
875  {
876  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
877  return -1;
878  }
879  int jj;
880  std::vector<int> vertList;
881  vertList.clear();
882  for (jj = 0; jj < int(snakein.isMeshVertIn.size()); ++jj)
883  {
884  if (snakein.isMeshVertIn[jj])
885  {
886  vertList.push_back(snakein.snakemesh()->verts(jj)->index);
887  }
888  }
889  if (int(snakein.isMeshVertIn.size()) == 0)
890  {
891  vertList.push_back(snakein.snakemesh()->verts(0)->index);
892  }
893  jj = PrintMesh(*(snakein.snakemesh()), strandID, timeStep, 4, vertList);
894  return (jj);
895 }
896 
897 // Functions triarray
898 void ExtractTriData(const triangulation &triout, triarray triangulation::*mp, int *nVert, int *nEdge, int *nVolu,
899  int *nSurf, int *totNumFaceNode)
900 {
901  // Extracts Data needed to write out a mesh to tecplot
902 
903  *nVert = (triout.*mp).size() * 3;
904  *nVolu = 1;
905  *nSurf = (triout.*mp).size();
906  *nEdge = (triout.*mp).size() * 3;
907  *totNumFaceNode = (triout.*mp).size() * 3;
908 }
909 
910 void ExtractTriData(int nTriangles, int *nVert, int *nEdge, int *nVolu, int *nSurf, int *totNumFaceNode)
911 {
912  // Extracts Data needed to write out a mesh to tecplot
913 
914  *nVert = nTriangles * 3;
915  *nVolu = 1;
916  *nSurf = nTriangles;
917  *nEdge = nTriangles * 3;
918  *totNumFaceNode = nTriangles * 3;
919 }
920 
921 int tecplotfile::VolDataBlock(const triangulation &triout, triarray triangulation::*mp, int nVert, int nVolu,
922  int nVertDat)
923 {
924  // Prints the Coord and Fill Data blocks to the tecplot file
925 
926  int ii, jj, kk;
927  int nTri, currInd, currType, nCoord;
928  nTri = (triout.*mp).size();
929  nCoord = int(triout.meshDep->verts(0)->coord.size());
930  // Print vertex Data
931  for (jj = 0; jj < nCoord; ++jj)
932  {
933  for (ii = 0; ii < nTri; ++ii)
934  {
935  for (kk = 0; kk < (3); ++kk)
936  {
937  currType = (triout.*mp)(ii)->pointtype[kk];
938  currInd = (triout.*mp)(ii)->pointind[kk];
939  if (currType == 1)
940  {
941  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
942  }
943  else if (currType == 2)
944  {
945  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
946  }
947  else if (currType == 3)
948  {
949  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
950  }
951  }
952  }
953  this->NewLine();
954  }
955  for (jj = nCoord; jj < nVertDat; ++jj)
956  {
957  for (ii = 0; ii < nVert; ++ii)
958  {
959  this->Print("%lf ", 0.0);
960  }
961  this->NewLine();
962  }
963 
964  // Print Cell Data
965  for (ii = 0; ii < nVolu; ++ii)
966  {
967  this->Print("%.16lf ", 0.0);
968  }
969  this->NewLine();
970  for (ii = 0; ii < nVolu; ++ii)
971  {
972  this->Print("%.16lf ", 0.0);
973  }
974  this->NewLine();
975  for (ii = 0; ii < nVolu; ++ii)
976  {
977  this->Print("%.16lf ", 0.0);
978  }
979  this->NewLine();
980  return (0);
981 }
982 
983 int tecplotfile::SurfDataBlock(const triangulation &triout, triarray triangulation::*mp, int nVert, int nSurf,
984  int nVertDat)
985 {
986  // Prints the Coord and Fill Data blocks to the tecplot file
987 
988  int ii, jj, kk;
989  int nTri, currInd, currType, nCoord;
990  nTri = (triout.*mp).size();
991  nCoord = int(triout.meshDep->verts(0)->coord.size());
992  // Print vertex Data
993  for (jj = 0; jj < nCoord; ++jj)
994  {
995  for (ii = 0; ii < nTri; ++ii)
996  {
997  for (kk = 0; kk < (3); ++kk)
998  {
999  currType = (triout.*mp)(ii)->pointtype[kk];
1000  currInd = (triout.*mp)(ii)->pointind[kk];
1001  if (currType == 1)
1002  {
1003  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
1004  }
1005  else if (currType == 2)
1006  {
1007  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
1008  }
1009  else if (currType == 3)
1010  {
1011  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
1012  }
1013  }
1014  }
1015  this->NewLine();
1016  }
1017  for (jj = nCoord; jj < nVertDat; ++jj)
1018  {
1019  for (ii = 0; ii < nVert; ++ii)
1020  {
1021  this->Print("%lf ", 0.0);
1022  }
1023  this->NewLine();
1024  }
1025 
1026  // Print Cell Data
1027  for (ii = 0; ii < nSurf; ++ii)
1028  {
1029  this->Print("%.16lf ", 0.0);
1030  }
1031  this->NewLine();
1032  for (ii = 0; ii < nSurf; ++ii)
1033  {
1034  this->Print("%.16lf ", 0.0);
1035  }
1036  this->NewLine();
1037  for (ii = 0; ii < nSurf; ++ii)
1038  {
1039  this->Print("%.16lf ", 0.0);
1040  }
1041  this->NewLine();
1042  return (0);
1043 }
1044 
1045 int tecplotfile::LineDataBlock(const triangulation &triout, triarray triangulation::*mp, int nVert, int nEdge,
1046  int nVertDat, int nCellDat)
1047 {
1048  // Prints the Coord and Fill Data blocks to the tecplot file
1049 
1050  int ii, jj, kk;
1051  int nTri, currInd, currType, nCoord;
1052  nTri = (triout.*mp).size();
1053  nCoord = int(triout.meshDep->verts(0)->coord.size());
1054  // Print vertex Data
1055  for (jj = 0; jj < nCoord; ++jj)
1056  {
1057  for (ii = 0; ii < nTri; ++ii)
1058  {
1059  for (kk = 0; kk < (3); ++kk)
1060  {
1061  currType = (triout.*mp)(ii)->pointtype[kk];
1062  currInd = (triout.*mp)(ii)->pointind[kk];
1063  if (currType == 1)
1064  {
1065  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
1066  }
1067  else if (currType == 2)
1068  {
1069  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
1070  }
1071  else if (currType == 3)
1072  {
1073  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
1074  }
1075  else
1076  {
1077  RSVS3D_ERROR_ARGUMENT("unknown point type");
1078  }
1079  }
1080  }
1081  this->NewLine();
1082  }
1083  for (jj = nCoord; jj < nVertDat; ++jj)
1084  {
1085  for (ii = 0; ii < nVert; ++ii)
1086  {
1087  this->Print("%lf ", 0.0);
1088  }
1089  this->NewLine();
1090  }
1091  // Print Cell Data
1092  for (jj = 0; jj < nCellDat; ++jj)
1093  {
1094  for (ii = 0; ii < nEdge; ++ii)
1095  {
1096  this->Print("%i ", 0);
1097  }
1098  this->NewLine();
1099  }
1100  return (0);
1101 }
1102 
1103 int tecplotfile::LineDataBlock(const triangulation &triout, triarray triangulation::*mp, int nVert, int nEdge,
1104  int nVertDat, int nCellDat, const vector<int> &triList)
1105 {
1106  // Prints the Coord and Fill Data blocks to the tecplot file
1107 
1108  int ii, jj, kk;
1109  int nTri, currInd, currType, nCoord;
1110  nTri = int(triList.size());
1111  nCoord = int(triout.meshDep->verts(0)->coord.size());
1112  // Print vertex Data
1113  for (jj = 0; jj < nCoord; ++jj)
1114  {
1115  for (ii = 0; ii < nTri; ++ii)
1116  {
1117  for (kk = 0; kk < (3); ++kk)
1118  {
1119  currType = (triout.*mp).isearch(triList[ii])->pointtype[kk];
1120  currInd = (triout.*mp).isearch(triList[ii])->pointind[kk];
1121  if (currType == 1)
1122  {
1123  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
1124  }
1125  else if (currType == 2)
1126  {
1127  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
1128  }
1129  else if (currType == 3)
1130  {
1131  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
1132  }
1133  else
1134  {
1135  RSVS3D_ERROR_ARGUMENT("unknown point type");
1136  }
1137  }
1138  }
1139  this->NewLine();
1140  }
1141  for (jj = nCoord; jj < nVertDat; ++jj)
1142  {
1143  for (ii = 0; ii < nVert; ++ii)
1144  {
1145  this->Print("%lf ", 0.0);
1146  }
1147  this->NewLine();
1148  }
1149  // Print Cell Data
1150  for (jj = 0; jj < nCellDat; ++jj)
1151  {
1152  for (ii = 0; ii < nEdge; ++ii)
1153  {
1154  this->Print("%i ", 0);
1155  }
1156  this->NewLine();
1157  }
1158  return (0);
1159 }
1160 
1161 int tecplotfile::SurfFaceMap(const triangulation &triout, triarray triangulation::*mp)
1162 {
1163  int ii, jj, kk;
1164 
1165  int nTri;
1166  nTri = (triout.*mp).size();
1167  // RSVS3D_ERROR_ARGUMENT("Surface Map not supported for triangulation")
1168  kk = 1;
1169  for (ii = 0; ii < nTri; ++ii)
1170  { // Print Number of vertices per face
1171  for (jj = 0; jj < 3; ++jj)
1172  {
1173  if (jj == 2)
1174  {
1175  this->Print("%i %i\n", kk - 2, kk);
1176  }
1177  else
1178  {
1179  this->Print("%i %i\n", kk, kk + 1);
1180  }
1181  this->ResetLine();
1182  ++kk;
1183  }
1184  }
1185 
1186  for (ii = 0; ii < nTri; ++ii)
1187  { // Print connectedFace number
1188  for (jj = 0; jj < 3; ++jj)
1189  {
1190  this->Print("%i ", ii + 1);
1191  }
1192  }
1193  this->NewLine();
1194  for (ii = 0; ii < 3 * nTri; ++ii)
1195  { // Print 0
1196  this->Print("%i ", 0);
1197  }
1198  this->NewLine();
1199  return (0);
1200 }
1201 
1202 int tecplotfile::LineFaceMap(const triangulation &triout, triarray triangulation::*mp)
1203 {
1204  int ii, jj, kk;
1205 
1206  int nTri;
1207  nTri = (triout.*mp).size();
1208  // RSVS3D_ERROR_ARGUMENT("Surface Map not supported for triangulation")
1209  kk = 1;
1210  for (ii = 0; ii < nTri; ++ii)
1211  { // Print Number of vertices per face
1212  for (jj = 0; jj < 3; ++jj)
1213  {
1214  if (jj == 2)
1215  {
1216  this->Print("%i %i\n", kk - 2, kk);
1217  }
1218  else
1219  {
1220  this->Print("%i %i\n", kk, kk + 1);
1221  }
1222  this->ResetLine();
1223  ++kk;
1224  }
1225  }
1226 
1227  return (0);
1228 }
1229 
1230 int tecplotfile::LineFaceMap(const vector<int> &triList)
1231 {
1232  int ii, jj, kk;
1233 
1234  int nTri;
1235  nTri = int(triList.size());
1236  // RSVS3D_ERROR_ARGUMENT("Surface Map not supported for triangulation")
1237  kk = 1;
1238  for (ii = 0; ii < nTri; ++ii)
1239  { // Print Number of vertices per face
1240  for (jj = 0; jj < 3; ++jj)
1241  {
1242  if (jj == 2)
1243  {
1244  this->Print("%i %i\n", kk - 2, kk);
1245  }
1246  else
1247  {
1248  this->Print("%i %i\n", kk, kk + 1);
1249  }
1250  this->ResetLine();
1251  ++kk;
1252  }
1253  }
1254 
1255  return (0);
1256 }
1257 
1258 int tecplotfile::VolFaceMap(const triangulation &triout, triarray triangulation::*mp, int nSurf)
1259 {
1260  int ii, jj, kk, n;
1261  n = (triout.*mp)(0)->pointind.size();
1262  for (ii = 0; ii < nSurf; ++ii)
1263  { // Print Number of vertices per face
1264  this->Print("%i ", n);
1265  }
1266  this->NewLine();
1267  kk = 1;
1268  for (ii = 0; ii < nSurf; ++ii)
1269  { // print ordered list of vertices in face
1270  for (jj = 0; jj < 3; ++jj)
1271  {
1272  this->Print("%i ", kk);
1273  kk++;
1274  }
1275  this->NewLine();
1276  }
1277  for (ii = 0; ii < nSurf; ++ii)
1278  {
1279  this->Print("%i ", 1);
1280  }
1281  this->NewLine();
1282  for (ii = 0; ii < nSurf; ++ii)
1283  {
1284  this->Print("%i ", 0);
1285  }
1286 
1287  return (0);
1288 }
1289 // Class function Implementation
1290 int tecplotfile::PrintTriangulation(const triangulation &triout, triarray triangulation::*mp, int strandID,
1291  double timeStep, int forceOutType, const vector<int> &triList)
1292 {
1293  if (!this->isOpen())
1294  {
1295  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
1296  return -1;
1297  }
1298  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat;
1299 
1300  if (triList.size() == 0 || forceOutType != 3)
1301  {
1302  ExtractTriData(triout, mp, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
1303  }
1304  else
1305  {
1306  ExtractTriData(triList.size(), &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
1307  }
1308  if (nVert > 0)
1309  {
1310  if (nZones == 0)
1311  {
1312  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"\n");
1313  }
1314 
1315  this->NewZone();
1316 
1317  if (strandID > 0)
1318  {
1319  this->StrandTime(strandID, timeStep);
1320  }
1321  // Fixed by the dimensionality of the mesh
1322  nVertDat = 3;
1323  nCellDat = 3;
1324 
1325  if (forceOutType == tecplotconst::autoselect)
1326  {
1327  if (nVolu > 0 && triout.snakeDep->Check3D())
1328  {
1329  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
1330  }
1331  else if (nSurf > 0)
1332  {
1333  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
1334  }
1335  else
1336  {
1337  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
1338  }
1339  }
1340 
1341  if (forceOutType == tecplotconst::polyhedron)
1342  {
1343  this->ZoneHeaderPolyhedron(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat);
1344  this->VolDataBlock(triout, mp, nVert, nVolu, nVertDat);
1345  this->VolFaceMap(triout, mp, nSurf);
1346  }
1347  else if (forceOutType == tecplotconst::polygon)
1348  {
1349  this->ZoneHeaderPolygon(nVert, nEdge, nSurf, nVertDat, nCellDat);
1350  this->SurfDataBlock(triout, mp, nVert, nSurf, nVertDat);
1351  this->SurfFaceMap(triout, mp);
1352  }
1353  else if (forceOutType == tecplotconst::line)
1354  {
1355  this->ZoneHeaderFelineseg(nVert, nEdge, nVertDat, nCellDat);
1356  if (triList.size() == 0)
1357  {
1358  this->LineDataBlock(triout, mp, nVert, nEdge, nVertDat, nCellDat);
1359  this->LineFaceMap(triout, mp);
1360  }
1361  else
1362  {
1363  this->LineDataBlock(triout, mp, nVert, nEdge, nVertDat, nCellDat, triList);
1364  this->LineFaceMap(triList);
1365  }
1366  }
1367  }
1368  return (0);
1369 }
1370 
1371 // Functions
1372 void ExtractTriData(const triangulation &triout, trisurfarray triangulation::*mp, int *nVert, int *nEdge, int *nVolu,
1373  int *nSurf, int *totNumFaceNode)
1374 {
1375  // Extracts Data needed to write out a mesh to tecplot
1376  int ii;
1377 
1378  *nSurf = (triout.*mp).size();
1379  *nVolu = 1;
1380  *nVert = 0;
1381  *nEdge = 0;
1382  for (ii = 0; ii < *nSurf; ii++)
1383  {
1384  *nVert += (triout.*mp)(ii)->indvert.size();
1385  }
1386 
1387  *nEdge = *nVert;
1388  *totNumFaceNode = *nVert;
1389 }
1390 
1391 int tecplotfile::VolDataBlock(const triangulation &triout, trisurfarray triangulation::*mp, int nVert, int nVolu,
1392  int nVertDat)
1393 {
1394  // Prints the Coord and Fill Data blocks to the tecplot file
1395 
1396  int ii, jj, kk;
1397  int nTri, currInd, currType, nCoord;
1398  nTri = (triout.*mp).size();
1399  nCoord = int(triout.meshDep->verts(0)->coord.size());
1400  // Print vertex Data
1401  for (jj = 0; jj < nCoord; ++jj)
1402  {
1403  for (ii = 0; ii < nTri; ++ii)
1404  {
1405  for (kk = 0; kk < int((triout.*mp)(ii)->typevert.size()); ++kk)
1406  {
1407  currType = (triout.*mp)(ii)->typevert[kk];
1408  currInd = (triout.*mp)(ii)->indvert[kk];
1409  if (currType == 1)
1410  {
1411  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
1412  }
1413  else if (currType == 2)
1414  {
1415  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
1416  }
1417  else if (currType == 3)
1418  {
1419  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
1420  }
1421  }
1422  }
1423  this->NewLine();
1424  }
1425  for (jj = nCoord; jj < nVertDat; ++jj)
1426  {
1427  for (ii = 0; ii < nVert; ++ii)
1428  {
1429  this->Print("%lf ", 0.0);
1430  }
1431  this->NewLine();
1432  }
1433 
1434  // Print Cell Data
1435  for (ii = 0; ii < nVolu; ++ii)
1436  {
1437  this->Print("%.16lf ", 0.0);
1438  }
1439  this->NewLine();
1440  for (ii = 0; ii < nVolu; ++ii)
1441  {
1442  this->Print("%.16lf ", 0.0);
1443  }
1444  this->NewLine();
1445  for (ii = 0; ii < nVolu; ++ii)
1446  {
1447  this->Print("%.16lf ", 0.0);
1448  }
1449  this->NewLine();
1450  return (0);
1451 }
1452 
1453 int tecplotfile::SurfDataBlock(const triangulation &triout, trisurfarray triangulation::*mp, int nVert, int nSurf,
1454  int nVertDat)
1455 {
1456  // Prints the Coord and Fill Data blocks to the tecplot file
1457 
1458  int ii, jj, kk;
1459  int nTri, currInd, currType, nCoord;
1460  nTri = (triout.*mp).size();
1461  nCoord = int(triout.meshDep->verts(0)->coord.size());
1462  // Print vertex Data
1463  for (jj = 0; jj < nCoord; ++jj)
1464  {
1465  for (ii = 0; ii < nTri; ++ii)
1466  {
1467  for (kk = 0; kk < int((triout.*mp)(ii)->typevert.size()); ++kk)
1468  {
1469  currType = (triout.*mp)(ii)->typevert[kk];
1470  currInd = (triout.*mp)(ii)->indvert[kk];
1471  if (currType == 1)
1472  {
1473  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
1474  }
1475  else if (currType == 2)
1476  {
1477  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
1478  }
1479  else if (currType == 3)
1480  {
1481  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
1482  }
1483  }
1484  }
1485  this->NewLine();
1486  }
1487  for (jj = nCoord; jj < nVertDat; ++jj)
1488  {
1489  for (ii = 0; ii < nVert; ++ii)
1490  {
1491  this->Print("%lf ", 0.0);
1492  }
1493  this->NewLine();
1494  }
1495 
1496  // Print Cell Data
1497  for (ii = 0; ii < nSurf; ++ii)
1498  {
1499  this->Print("%.16lf ", 0.0);
1500  }
1501  this->NewLine();
1502  for (ii = 0; ii < nSurf; ++ii)
1503  {
1504  this->Print("%.16lf ", 0.0);
1505  }
1506  this->NewLine();
1507  for (ii = 0; ii < nSurf; ++ii)
1508  {
1509  this->Print("%.16lf ", 0.0);
1510  }
1511  this->NewLine();
1512  return (0);
1513 }
1514 
1515 int tecplotfile::LineDataBlock(const triangulation &triout, trisurfarray triangulation::*mp, int nVert, int nEdge,
1516  int nVertDat, int nCellDat)
1517 {
1518  // Prints the Coord and Fill Data blocks to the tecplot file
1519 
1520  int ii, jj, kk;
1521  int nTri, currInd, currType, nCoord;
1522  nTri = (triout.*mp).size();
1523  nCoord = int(triout.meshDep->verts(0)->coord.size());
1524  // Print vertex Data
1525  for (jj = 0; jj < nCoord; ++jj)
1526  {
1527  for (ii = 0; ii < nTri; ++ii)
1528  {
1529  for (kk = 0; kk < int((triout.*mp)(ii)->typevert.size()); ++kk)
1530  {
1531  currType = (triout.*mp)(ii)->typevert[kk];
1532  currInd = (triout.*mp)(ii)->indvert[kk];
1533  if (currType == 1)
1534  {
1535  this->Print("%.16lf ", triout.meshDep->verts.isearch(currInd)->coord[jj]);
1536  }
1537  else if (currType == 2)
1538  {
1539  this->Print("%.16lf ", triout.snakeDep->snakeconn.verts.isearch(currInd)->coord[jj]);
1540  }
1541  else if (currType == 3)
1542  {
1543  this->Print("%.16lf ", triout.trivert.isearch(currInd)->coord(jj));
1544  }
1545  else
1546  {
1547  RSVS3D_ERROR_ARGUMENT("unknown point type");
1548  }
1549  }
1550  }
1551  this->NewLine();
1552  }
1553  for (jj = nCoord; jj < nVertDat; ++jj)
1554  {
1555  for (ii = 0; ii < nVert; ++ii)
1556  {
1557  this->Print("%lf ", 0.0);
1558  }
1559  this->NewLine();
1560  }
1561  // Print Cell Data
1562  for (jj = 0; jj < nCellDat; ++jj)
1563  {
1564  for (ii = 0; ii < nEdge; ++ii)
1565  {
1566  this->Print("%i ", 0);
1567  }
1568  this->NewLine();
1569  }
1570  return (0);
1571 }
1572 int tecplotfile::SurfFaceMap(const triangulation &triout, trisurfarray triangulation::*mp)
1573 {
1574  int ii, jj, kk;
1575 
1576  int nTri;
1577  nTri = (triout.*mp).size();
1578  // RSVS3D_ERROR_ARGUMENT("Surface Map not supported for triangulation")
1579  kk = 1;
1580  for (ii = 0; ii < nTri; ++ii)
1581  { // Print Number of vertices per face
1582  for (jj = 0; jj < int((triout.*mp)(ii)->typevert.size()); ++jj)
1583  {
1584  if (jj == (int((triout.*mp)(ii)->typevert.size()) - 1))
1585  {
1586  this->Print("%i %i\n", kk - (int((triout.*mp)(ii)->typevert.size()) - 1), kk);
1587  }
1588  else
1589  {
1590  this->Print("%i %i\n", kk, kk + 1);
1591  }
1592  this->ResetLine();
1593  ++kk;
1594  }
1595  }
1596 
1597  for (ii = 0; ii < nTri; ++ii)
1598  { // Print connectedFace number
1599  for (jj = 0; jj < 3; ++jj)
1600  {
1601  this->Print("%i ", ii + 1);
1602  }
1603  }
1604  this->NewLine();
1605  for (ii = 0; ii < 3 * nTri; ++ii)
1606  { // Print 0
1607  this->Print("%i ", 0);
1608  }
1609  this->NewLine();
1610  return (0);
1611 }
1612 
1613 int tecplotfile::LineFaceMap(const triangulation &triout, trisurfarray triangulation::*mp)
1614 {
1615  int ii, jj, kk;
1616 
1617  int nTri;
1618  nTri = (triout.*mp).size();
1619  // RSVS3D_ERROR_ARGUMENT("Surface Map not supported for triangulation")
1620  kk = 1;
1621  for (ii = 0; ii < nTri; ++ii)
1622  { // Print Number of vertices per face
1623  for (jj = 0; jj < int((triout.*mp)(ii)->typevert.size()); ++jj)
1624  {
1625  if (jj == (int((triout.*mp)(ii)->typevert.size()) - 1))
1626  {
1627  this->Print("%i %i\n", kk - (int((triout.*mp)(ii)->typevert.size()) - 1), kk);
1628  }
1629  else
1630  {
1631  this->Print("%i %i\n", kk, kk + 1);
1632  }
1633  this->ResetLine();
1634  ++kk;
1635  }
1636  }
1637 
1638  return (0);
1639 }
1640 
1641 int tecplotfile::VolFaceMap(const triangulation &triout, trisurfarray triangulation::*mp, int nSurf)
1642 {
1643  int ii, jj, kk;
1644 
1645  for (ii = 0; ii < nSurf; ++ii)
1646  { // Print Number of vertices per face
1647  this->Print("%i ", int((triout.*mp)(ii)->typevert.size()));
1648  }
1649  this->NewLine();
1650  kk = 1;
1651  for (ii = 0; ii < nSurf; ++ii)
1652  { // print ordered list of vertices in face
1653  for (jj = 0; jj < int((triout.*mp)(ii)->typevert.size()); ++jj)
1654  {
1655  this->Print("%i ", kk);
1656  kk++;
1657  }
1658  this->NewLine();
1659  }
1660  for (ii = 0; ii < nSurf; ++ii)
1661  {
1662  this->Print("%i ", 1);
1663  }
1664  this->NewLine();
1665  for (ii = 0; ii < nSurf; ++ii)
1666  {
1667  this->Print("%i ", 0);
1668  }
1669 
1670  return (0);
1671 }
1672 // Class function Implementation
1673 int tecplotfile::PrintTriangulation(const triangulation &triout, trisurfarray triangulation::*mp, int strandID,
1674  double timeStep, int forceOutType)
1675 {
1676  if (!this->isOpen())
1677  {
1678  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
1679  return -1;
1680  }
1681  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat;
1682 
1683  ExtractTriData(triout, mp, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
1684  if (nVert > 0)
1685  {
1686  if (nZones == 0)
1687  {
1688  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"\n");
1689  }
1690 
1691  this->NewZone();
1692 
1693  if (strandID > 0)
1694  {
1695  this->StrandTime(strandID, timeStep);
1696  }
1697  // Fixed by the dimensionality of the mesh
1698  nVertDat = 3;
1699  nCellDat = 3;
1700 
1701  if (forceOutType == tecplotconst::autoselect)
1702  {
1703  if (nVolu > 0 && triout.snakeDep->Check3D())
1704  {
1705  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
1706  }
1707  else if (nSurf > 0)
1708  {
1709  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
1710  }
1711  else
1712  {
1713  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
1714  }
1715  }
1716 
1717  if (forceOutType == tecplotconst::polyhedron)
1718  {
1719  this->ZoneHeaderPolyhedron(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat);
1720  this->VolDataBlock(triout, mp, nVert, nVolu, nVertDat);
1721  this->VolFaceMap(triout, mp, nSurf);
1722  }
1723  else if (forceOutType == tecplotconst::polygon)
1724  {
1725  this->ZoneHeaderPolygon(nVert, nEdge, nSurf, nVertDat, nCellDat);
1726  this->SurfDataBlock(triout, mp, nVert, nSurf, nVertDat);
1727  this->SurfFaceMap(triout, mp);
1728  }
1729  else if (forceOutType == tecplotconst::line)
1730  {
1731  this->ZoneHeaderFelineseg(nVert, nEdge, nVertDat, nCellDat);
1732  this->LineDataBlock(triout, mp, nVert, nEdge, nVertDat, nCellDat);
1733  this->LineFaceMap(triout, mp);
1734  }
1735  }
1736  return (0);
1737 }
1738 
1739 void tecplotfile::ZoneHeaderPolyhedron(int nVert, int nVolu, int nSurf, int totNumFaceNode, int nVertDat, int nCellDat)
1740 {
1741  this->Print("ZONETYPE = FEPOLYHEDRON\n");
1742  this->Print("NODES = %i\n", nVert);
1743  this->Print("ELEMENTS = %i\n", nVolu);
1744  this->Print("FACES = %i\n", nSurf);
1745  this->Print("TOTALNUMFACENODES = %i\n", totNumFaceNode);
1746  this->Print("NUMCONNECTEDBOUNDARYFACES = 0\n");
1747  this->Print("TOTALNUMBOUNDARYCONNECTIONS = 0\n");
1748  this->Print("VARLOCATION=([%i-%i]=NODAL ,[%i-%i]=CELLCENTERED)\n", 1, nVertDat, nVertDat + 1, nVertDat + nCellDat);
1749  this->Print("DATAPACKING=BLOCK\n");
1750 }
1751 
1752 void tecplotfile::ZoneHeaderPolygon(int nVert, int nEdge, int nSurf, int nVertDat, int nCellDat)
1753 {
1754  this->Print("ZONETYPE = FEPOLYGON\n");
1755  this->Print("NODES = %i\n", nVert);
1756  this->Print("ELEMENTS = %i\n", nSurf);
1757  this->Print("FACES = %i\n", nEdge);
1758  this->Print("TOTALNUMFACENODES = %i\n", 2 * nEdge);
1759  this->Print("NUMCONNECTEDBOUNDARYFACES = 0\n");
1760  this->Print("TOTALNUMBOUNDARYCONNECTIONS = 0\n");
1761  this->Print("VARLOCATION=([%i-%i]=NODAL ,[%i-%i]=CELLCENTERED)\n", 1, nVertDat, nVertDat + 1, nVertDat + nCellDat);
1762  this->Print("DATAPACKING=BLOCK\n");
1763 }
1764 
1765 void tecplotfile::ZoneHeaderFelineseg(int nVert, int nEdge, int nVertDat, int nCellDat)
1766 {
1767  this->Print("ZONETYPE = FELINESEG\n");
1768  this->Print("NODES = %i\n", nVert);
1769  this->Print("ELEMENTS = %i\n", nEdge);
1770  this->Print("FACES = %i\n", nEdge);
1771  this->Print("VARLOCATION=([%i-%i]=NODAL ,[%i-%i]=CELLCENTERED)\n", 1, nVertDat, nVertDat + 1, nVertDat + nCellDat);
1772  this->Print("DATAPACKING=BLOCK\n");
1773 }
1774 
1775 void tecplotfile::ZoneHeaderPolyhedronSnake(int nVert, int nVolu, int nSurf, int totNumFaceNode, int nVertDat,
1776  int nCellDat, int nSensDat)
1777 {
1778  this->Print("ZONETYPE = FEPOLYHEDRON\n");
1779  this->Print("NODES = %i\n", nVert);
1780  this->Print("ELEMENTS = %i\n", nVolu);
1781  this->Print("FACES = %i\n", nSurf);
1782  this->Print("TOTALNUMFACENODES = %i\n", totNumFaceNode);
1783  this->Print("NUMCONNECTEDBOUNDARYFACES = 0\n");
1784  this->Print("TOTALNUMBOUNDARYCONNECTIONS = 0\n");
1785  this->Print("VARLOCATION=([%i-%i]=NODAL)\n", 1, nVertDat + nCellDat + nSensDat);
1786 
1787  this->Print("DATAPACKING=BLOCK\n");
1788 }
1789 
1790 void tecplotfile::ZoneHeaderPolygonSnake(int nVert, int nEdge, int nSurf, int nVertDat, int nCellDat, int nSensDat)
1791 {
1792  this->Print("ZONETYPE = FEPOLYGON\n");
1793  this->Print("NODES = %i\n", nVert);
1794  this->Print("ELEMENTS = %i\n", nSurf);
1795  this->Print("FACES = %i\n", nEdge);
1796  this->Print("TOTALNUMFACENODES = %i\n", 2 * nEdge);
1797  this->Print("NUMCONNECTEDBOUNDARYFACES = 0\n");
1798  this->Print("TOTALNUMBOUNDARYCONNECTIONS = 0\n");
1799  this->Print("VARLOCATION=([%i-%i]=NODAL)\n", 1, nVertDat + nCellDat + nSensDat);
1800 
1801  this->Print("DATAPACKING=BLOCK\n");
1802 }
1803 
1804 void tecplotfile::ZoneHeaderFelinesegSnake(int nVert, int nEdge, int nVertDat, int nCellDat, int nSensDat)
1805 {
1806  this->Print("ZONETYPE = FELINESEG\n");
1807  this->Print("NODES = %i\n", nVert);
1808  this->Print("ELEMENTS = %i\n", nEdge);
1809  this->Print("FACES = %i\n", nEdge);
1810  this->Print("VARLOCATION=([%i-%i]=NODAL)\n", 1, nVertDat + nCellDat + nSensDat);
1811 
1812  this->Print("DATAPACKING=BLOCK\n");
1813 }
1814 
1815 void tecplotfile::ZoneHeaderOrdered(int nVert, int nVertDat, int nCellDat, int nSensDat)
1816 {
1817  this->Print("ZONETYPE = ORDERED\n");
1818  this->Print("I = %i\n", nVert);
1819  this->Print("VARLOCATION=([%i-%i]=NODAL)\n", 1, nVertDat + nCellDat + nSensDat);
1820  this->Print("DATAPACKING=BLOCK\n");
1821 }
1822 
1823 // Print snake with sensitivity
1824 int tecplotfile::PrintSnakeSensitivity(const triangulation &triRSVS, const RSVScalc &calcObj, int strandID,
1825  double timeStep, int forceOutType, const vector<int> &vertList)
1826 {
1827  if (!this->isOpen())
1828  {
1829  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
1830  return -1;
1831  }
1832  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat;
1833  const snake &snakeout = *triRSVS.snakeDep;
1834  ExtractMeshData(snakeout.snakeconn, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
1835  if (nVert == 0)
1836  { // Don't print mesh
1837  return (1);
1838  }
1839 
1840  if (nZones == 0)
1841  {
1842  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"");
1843  for (int i = 0; i < calcObj.numConstr(); ++i)
1844  {
1845  this->Print(", \"sens_%i\"", i);
1846  }
1847  this->Print("\n");
1848  }
1849 
1850  this->NewZone();
1851 
1852  if (strandID > 0)
1853  {
1854  this->StrandTime(strandID, timeStep);
1855  }
1856  // Fixed by the dimensionality of the mesh
1857  nVertDat = 3;
1858  nCellDat = 3;
1859  nSensDat = calcObj.numConstr();
1860 
1861  if (forceOutType == tecplotconst::autoselect)
1862  {
1863  if (nVolu > 0)
1864  {
1865  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
1866  }
1867  else if (nSurf > 0)
1868  {
1869  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
1870  }
1871  else if (nEdge > 0)
1872  {
1873  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
1874  }
1875  else
1876  {
1877  forceOutType = tecplotconst::point;
1878  }
1879  }
1880 
1881  if (forceOutType == tecplotconst::polyhedron)
1882  {
1883  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat);
1884  this->SnakeDataBlock(snakeout, nVert, nVertDat);
1885  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
1886  this->VolFaceMap(snakeout.snakeconn, nSurf);
1887  }
1888  else if (forceOutType == tecplotconst::polygon)
1889  {
1890  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat, nSensDat);
1891  this->SnakeDataBlock(snakeout, nVert, nVertDat);
1892  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
1893  this->SurfFaceMap(snakeout.snakeconn, nEdge);
1894  }
1895  else if (forceOutType == tecplotconst::line)
1896  {
1897  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat, nSensDat);
1898  this->SnakeDataBlock(snakeout, nVert, nVertDat);
1899  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
1900  this->LineFaceMap(snakeout.snakeconn, nEdge);
1901  }
1902  else if (forceOutType == tecplotconst::point)
1903  {
1904  if (int(vertList.size()) == nVert)
1905  {
1906  nVert = 0;
1907  for (int ii = 0; ii < int(vertList.size()); ++ii)
1908  {
1909  nVert += int(vertList[ii]);
1910  }
1911  }
1912  else if (vertList.size() > 0)
1913  {
1914  nVert = int(vertList.size());
1915  }
1916  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat, nSensDat);
1917  this->SnakeDataBlock(snakeout, nVert, nVertDat);
1918  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
1919  // No map just points
1920  }
1921  return (0);
1922 }
1923 
1924 int tecplotfile::PrintSnakeGradients(const triangulation &triRSVS, const RSVScalc &calcObj, int strandID,
1925  double timeStep, int forceOutType, const vector<int> &vertList)
1926 {
1927  if (!this->isOpen())
1928  {
1929  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
1930  return -1;
1931  }
1932  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat;
1933  const snake &snakeout = *triRSVS.snakeDep;
1934  ExtractMeshData(snakeout.snakeconn, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
1935  if (nVert == 0)
1936  { // Don't print mesh
1937  return (1);
1938  }
1939  int nPreConstr = 6;
1940  if (nZones == 0)
1941  {
1942  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"");
1943  this->Print(", \"dObj\"");
1944  this->Print(", \"dLag\"");
1945  this->Print(", \"HObj\"");
1946  this->Print(", \"HConstr\"");
1947  this->Print(", \"HLag\"");
1948  this->Print(", \"deltaDV\"");
1949  for (int i = 0; i < calcObj.numConstr(); ++i)
1950  {
1951  this->Print(", \"dConstr_%i\"", i);
1952  }
1953  this->Print("\n");
1954  }
1955 
1956  this->NewZone();
1957 
1958  if (strandID > 0)
1959  {
1960  this->StrandTime(strandID, timeStep);
1961  }
1962  // Fixed by the dimensionality of the mesh
1963  nVertDat = 3;
1964  nCellDat = 3;
1965  nSensDat = calcObj.numConstr() + nPreConstr;
1966 
1967  if (forceOutType == tecplotconst::autoselect)
1968  {
1969  if (nVolu > 0)
1970  {
1971  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
1972  }
1973  else if (nSurf > 0)
1974  {
1975  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
1976  }
1977  else if (nEdge > 0)
1978  {
1979  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
1980  }
1981  else
1982  {
1983  forceOutType = tecplotconst::point;
1984  }
1985  }
1986 
1987  if (forceOutType == tecplotconst::polyhedron)
1988  {
1989  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat);
1990  this->SnakeDataBlock(snakeout, nVert, nVertDat);
1991  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat - nPreConstr, -nPreConstr, 2);
1992  this->VolFaceMap(snakeout.snakeconn, nSurf);
1993  }
1994  else if (forceOutType == tecplotconst::polygon)
1995  {
1996  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat, nSensDat);
1997  this->SnakeDataBlock(snakeout, nVert, nVertDat);
1998  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat - nPreConstr, -nPreConstr, 2);
1999  this->SurfFaceMap(snakeout.snakeconn, nEdge);
2000  }
2001  else if (forceOutType == tecplotconst::line)
2002  {
2003  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat, nSensDat);
2004  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2005  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat - nPreConstr, -nPreConstr, 2);
2006  this->LineFaceMap(snakeout.snakeconn, nEdge);
2007  }
2008  else if (forceOutType == tecplotconst::point)
2009  {
2010  if (int(vertList.size()) == nVert)
2011  {
2012  nVert = 0;
2013  for (int ii = 0; ii < int(vertList.size()); ++ii)
2014  {
2015  nVert += int(vertList[ii]);
2016  }
2017  }
2018  else if (vertList.size() > 0)
2019  {
2020  nVert = int(vertList.size());
2021  }
2022  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat, nSensDat);
2023  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2024  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat - nPreConstr, -nPreConstr, 2);
2025  // No map just points
2026  }
2027  return (0);
2028 }
2029 
2030 int tecplotfile::RSVScalcDataBlock(const triangulation &triRSVS, const RSVScalc &calcObj, int nVert, int nSensDat,
2031  int sensStart, int methodProcess)
2032 {
2033  if (!this->isOpen())
2034  {
2035  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
2036  return -1;
2037  }
2038  // Prints the Coord and Fill Data blocks to the tecplot file
2039 
2040  int ii, jj;
2041  // Print vertex Data
2042  std::vector<double> sensTemp;
2043  void (RSVScalc::*mp)(const triangulation &, std::vector<double> &, int) const;
2044  if (methodProcess == 1)
2045  {
2046  mp = &RSVScalc::ReturnSensitivities;
2047  }
2048  else if (methodProcess == 2)
2049  {
2050  mp = &RSVScalc::ReturnGradient;
2051  }
2052  else
2053  {
2054  RSVS3D_ERROR_ARGUMENT("Unknown methodProcess value accepted are 1 "
2055  "(sensitivities) and 2 (derivatives).");
2056  }
2057 
2058  for (jj = sensStart; jj < nSensDat; ++jj)
2059  {
2060  (calcObj.*mp)(triRSVS, sensTemp, jj);
2061  for (ii = 0; ii < nVert; ++ii)
2062  {
2063  this->Print("%.16lf ", sensTemp[ii]);
2064  }
2065  this->NewLine();
2066  }
2067  return 0;
2068 }
2069 
2070 int tecplotfile::RSVScalcVectorDataBlock(const triangulation &triRSVS, const RSVScalc &calcObj, int nVert,
2071  int numConstrPlot, int methodProcess)
2072 {
2073  if (!this->isOpen())
2074  {
2075  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
2076  return -1;
2077  }
2078  // Prints the Coord and Fill Data blocks to the tecplot file
2079 
2080  int ii;
2081  // Print vertex Data
2082  std::vector<double> sensTemp;
2083 
2084  void (RSVScalc::*mp)(const triangulation &, std::vector<double> &, int) const;
2085  if (methodProcess == 1)
2086  {
2087  mp = &RSVScalc::ReturnSensitivities;
2088  }
2089  else if (methodProcess == 2)
2090  {
2091  mp = &RSVScalc::ReturnGradient;
2092  }
2093  else
2094  {
2095  RSVS3D_ERROR_ARGUMENT("Unknown methodProcess value accepted are 1 "
2096  "(sensitivities) and 2 (derivatives).");
2097  }
2098  auto normals = MeshUnitNormals(triRSVS.snakeDep->snakeconn);
2099  auto directions = triRSVS.snakeDep->MoveDirections();
2100 
2101  (calcObj.*mp)(triRSVS, sensTemp, numConstrPlot);
2102  for (ii = 0; ii < nVert; ++ii)
2103  {
2104  double dotProd = 0.0;
2105  for (int j = 0; j < 3; ++j)
2106  {
2107  dotProd += normals[ii * 3 + j] * directions[ii * 3 + j];
2108  }
2109  this->Print("%.16lf ", sensTemp[ii]);
2110  sensTemp[ii] *= dotProd;
2111  }
2112  this->NewLine();
2113  for (int j = 0; j < 3; ++j)
2114  {
2115  for (ii = 0; ii < nVert; ++ii)
2116  {
2117  this->Print("%.16lf ", sensTemp[ii] * normals[ii * 3 + j]);
2118  }
2119  this->NewLine();
2120  }
2121 
2122  return 0;
2123 }
2124 
2125 // Print snake with sensitivity
2126 int tecplotfile::PrintSnakeSensitivityTime(const triangulation &triRSVS, const RSVScalc &calcObj, int strandID,
2127  double timeStep, int forceOutType, const vector<int> &vertList)
2128 {
2129  if (!this->isOpen())
2130  {
2131  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
2132  return -1;
2133  }
2134  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat, nSensStep;
2135  const snake &snakeout = *triRSVS.snakeDep;
2136  ExtractMeshData(snakeout.snakeconn, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
2137  if (nVert == 0)
2138  { // Don't print mesh
2139  return (1);
2140  }
2141  nVertDat = 3;
2142  nCellDat = 3;
2143  nSensDat = 1;
2144  nSensStep = calcObj.numConstr();
2145 
2146  if (nZones == 0)
2147  {
2148  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"");
2149  for (int i = 0; i < nSensDat; ++i)
2150  {
2151  this->Print(", \"sens_%i\"", i);
2152  }
2153  this->Print("\n");
2154  }
2155 
2156  this->NewZone();
2157  if (strandID <= 0)
2158  {
2159  while (this->isStrand(++strandID) && strandID <= 100)
2160  {
2161  }
2162  }
2163  if (strandID > 0)
2164  {
2165  this->StrandTime(strandID, timeStep);
2166  }
2167  // Fixed by the dimensionality of the mesh
2168 
2169  if (forceOutType == tecplotconst::autoselect)
2170  {
2171  if (nVolu > 0)
2172  {
2173  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
2174  }
2175  else if (nSurf > 0)
2176  {
2177  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
2178  }
2179  else if (nEdge > 0)
2180  {
2181  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
2182  }
2183  else
2184  {
2185  forceOutType = tecplotconst::point;
2186  }
2187  }
2188 
2189  if (forceOutType == tecplotconst::polyhedron)
2190  {
2191  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat);
2192  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2193  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
2194  this->VolFaceMap(snakeout.snakeconn, nSurf);
2195  }
2196  else if (forceOutType == tecplotconst::polygon)
2197  {
2198  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat, nSensDat);
2199  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2200  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
2201  this->SurfFaceMap(snakeout.snakeconn, nEdge);
2202  }
2203  else if (forceOutType == tecplotconst::line)
2204  {
2205  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat, nSensDat);
2206  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2207  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
2208  this->LineFaceMap(snakeout.snakeconn, nEdge);
2209  }
2210  else if (forceOutType == tecplotconst::point)
2211  {
2212  if (int(vertList.size()) == nVert)
2213  {
2214  nVert = 0;
2215  for (int ii = 0; ii < int(vertList.size()); ++ii)
2216  {
2217  nVert += int(vertList[ii]);
2218  }
2219  }
2220  else if (vertList.size() > 0)
2221  {
2222  nVert = int(vertList.size());
2223  }
2224  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat, nSensDat);
2225  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2226  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, nSensDat);
2227  // No map just points
2228  }
2229  int shareZone = this->ZoneNum();
2230  double tStepMultiplier = pow(10, -ceil(log10(double(nSensStep))));
2231  for (int i = nSensDat; i < nSensStep; ++i)
2232  {
2233  this->NewZone();
2234  this->StrandTime(strandID, timeStep + tStepMultiplier * i);
2235  if (forceOutType == tecplotconst::polyhedron)
2236  {
2237  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat);
2238  }
2239  else if (forceOutType == tecplotconst::polygon)
2240  {
2241  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat, nSensDat);
2242  }
2243  else if (forceOutType == tecplotconst::line)
2244  {
2245  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat, nSensDat);
2246  }
2247  else if (forceOutType == tecplotconst::point)
2248  {
2249  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat, nSensDat);
2250  }
2251  this->DefShareZoneVolume(shareZone, nVertDat + nCellDat);
2252  this->RSVScalcDataBlock(triRSVS, calcObj, nVert, i + 1, i);
2253  }
2254 
2255  return (0);
2256 }
2257 
2258 // Print snake with sensitivity
2259 int tecplotfile::PrintSnakeSensitivityVector(const triangulation &triRSVS, const RSVScalc &calcObj, int strandID,
2260  double timeStep, int forceOutType, const vector<int> &vertList)
2261 {
2262  if (!this->isOpen())
2263  {
2264  RSVS3D_ERROR_NOTHROW("File not open, no data printed.");
2265  return -1;
2266  }
2267  int nVert, nEdge, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat, nSensStep;
2268  const snake &snakeout = *triRSVS.snakeDep;
2269  ExtractMeshData(snakeout.snakeconn, &nVert, &nEdge, &nVolu, &nSurf, &totNumFaceNode);
2270  if (nVert == 0)
2271  { // Don't print mesh
2272  return (1);
2273  }
2274  nVertDat = 3;
2275  nCellDat = 3;
2276  nSensDat = 4;
2277  nSensStep = calcObj.numConstr();
2278 
2279  if (nZones == 0)
2280  {
2281  this->Print("VARIABLES = \"X\" ,\"Y\" , \"Z\" ,\"v1\" ,\"v2\", \"v3\"");
2282  for (int i = 0; i < nSensDat; ++i)
2283  {
2284  this->Print(", \"sens_%i\"", i);
2285  }
2286  this->Print("\n");
2287  }
2288 
2289  this->NewZone();
2290  if (strandID <= 0)
2291  {
2292  while (this->isStrand(++strandID) && strandID <= 100)
2293  {
2294  }
2295  }
2296  if (strandID > 0)
2297  {
2298  this->StrandTime(strandID, timeStep);
2299  }
2300  // Fixed by the dimensionality of the mesh
2301 
2302  if (forceOutType == tecplotconst::autoselect)
2303  {
2304  if (nVolu > 0)
2305  {
2306  forceOutType = tecplotconst::polyhedron; // output as volume data (FEPOLYHEDRON)
2307  }
2308  else if (nSurf > 0)
2309  {
2310  forceOutType = tecplotconst::polygon; // output as Surface data (FEPOLYGON)
2311  }
2312  else if (nEdge > 0)
2313  {
2314  forceOutType = tecplotconst::line; // output as line data (FELINESEG)
2315  }
2316  else
2317  {
2318  forceOutType = tecplotconst::point;
2319  }
2320  }
2321 
2322  if (forceOutType == tecplotconst::polyhedron)
2323  {
2324  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat);
2325  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2326  this->RSVScalcVectorDataBlock(triRSVS, calcObj, nVert);
2327  this->VolFaceMap(snakeout.snakeconn, nSurf);
2328  }
2329  else if (forceOutType == tecplotconst::polygon)
2330  {
2331  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat, nSensDat);
2332  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2333  this->RSVScalcVectorDataBlock(triRSVS, calcObj, nVert);
2334  this->SurfFaceMap(snakeout.snakeconn, nEdge);
2335  }
2336  else if (forceOutType == tecplotconst::line)
2337  {
2338  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat, nSensDat);
2339  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2340  this->RSVScalcVectorDataBlock(triRSVS, calcObj, nVert);
2341  this->LineFaceMap(snakeout.snakeconn, nEdge);
2342  }
2343  else if (forceOutType == tecplotconst::point)
2344  {
2345  if (int(vertList.size()) == nVert)
2346  {
2347  nVert = 0;
2348  for (int ii = 0; ii < int(vertList.size()); ++ii)
2349  {
2350  nVert += int(vertList[ii]);
2351  }
2352  }
2353  else if (vertList.size() > 0)
2354  {
2355  nVert = int(vertList.size());
2356  }
2357  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat, nSensDat);
2358  this->SnakeDataBlock(snakeout, nVert, nVertDat);
2359  this->RSVScalcVectorDataBlock(triRSVS, calcObj, nVert);
2360  // No map just points
2361  }
2362  int shareZone = this->ZoneNum();
2363  double tStepMultiplier = pow(10, -ceil(log10(double(nSensStep))));
2364  for (int i = 1; i < nSensStep; ++i)
2365  {
2366  this->NewZone();
2367  this->StrandTime(strandID, timeStep + tStepMultiplier * i);
2368  if (forceOutType == tecplotconst::polyhedron)
2369  {
2370  this->ZoneHeaderPolyhedronSnake(nVert, nVolu, nSurf, totNumFaceNode, nVertDat, nCellDat, nSensDat);
2371  }
2372  else if (forceOutType == tecplotconst::polygon)
2373  {
2374  this->ZoneHeaderPolygonSnake(nVert, nEdge, nSurf, nVertDat, nCellDat, nSensDat);
2375  }
2376  else if (forceOutType == tecplotconst::line)
2377  {
2378  this->ZoneHeaderFelinesegSnake(nVert, nEdge, nVertDat, nCellDat, nSensDat);
2379  }
2380  else if (forceOutType == tecplotconst::point)
2381  {
2382  this->ZoneHeaderOrdered(nVert, nVertDat, nCellDat, nSensDat);
2383  }
2384  this->DefShareZoneVolume(shareZone, nVertDat + nCellDat);
2385  this->RSVScalcVectorDataBlock(triRSVS, calcObj, nVert, i);
2386  }
2387 
2388  return (0);
2389 }
2390 
2391 // tecplotfile operations
2392 
2393 int tecplotfile::OpenFile(const char *str, const char *mode)
2394 {
2395  if (fid != NULL)
2396  {
2397  fclose(fid);
2398  }
2399  fid = fopen(str, mode);
2400  if (fid == NULL)
2401  {
2402  cout << "File '" << str << "' failed to open" << endl;
2403  return (-1);
2404  }
2405  return (0);
2406 }
2407 void tecplotfile::CloseFile()
2408 {
2409  if (fid != NULL)
2410  {
2411  fclose(fid);
2412  fid = NULL;
2413  }
2414 }
2415 
2416 // Test Code
2417 int Test_tecplotfile()
2418 {
2419  const char *fileToOpen;
2420  tecplotfile outmesh;
2421  int errFlag;
2422 
2423  fileToOpen = "../TESTOUT/tecout.plt";
2424 
2425  errFlag = outmesh.OpenFile(fileToOpen);
2426  if (errFlag != 0)
2427  {
2428  return (errFlag);
2429  }
2430  return (0);
2431 }
2432 
2433 int TestCompareReadWrite(const char *fileToOpen, mesh &blockGrid, tecplotfile &outmesh1)
2434 {
2435  int errFlag = 0;
2436  int errTest = 0;
2437  mesh blockGrid2;
2438  FILE *fid;
2439 
2440  fid = fopen(fileToOpen, "w");
2441  if (fid != NULL)
2442  {
2443  blockGrid.PrepareForUse();
2444  blockGrid.SetBorders();
2445  blockGrid.write(fid);
2446  fclose(fid);
2447  fid = fopen(fileToOpen, "r");
2448  if (fid != NULL)
2449  {
2450  blockGrid2.read(fid);
2451  fclose(fid);
2452  blockGrid.PrepareForUse();
2453  blockGrid2.PrepareForUse();
2454 
2455  outmesh1.PrintMesh(blockGrid2);
2456  errTest = CompareDisp(blockGrid, blockGrid2);
2457  if (!errTest)
2458  {
2459  // blockGrid.disp();
2460  // blockGrid2.disp();
2461  cerr << "Error: Displays were not the same after write read" << endl;
2462  errFlag++;
2463  }
2464  }
2465  else
2466  {
2467  cout << "File for mesh out failed to open" << endl;
2468  }
2469  }
2470  else
2471  {
2472  cout << "File for mesh out failed to open" << endl;
2473  }
2474  return (errFlag);
2475 }
RSVScalc.hpp
Provides the infrastructure for calculation of the RSVS equations.
RSVScalc
Class to handle the RSVS calculation.
Definition: RSVScalc.hpp:130
RSVScalc::numConstr
int numConstr() const
Getter for the number of constraints.
Definition: RSVScalc.hpp:455
mesh
Class for mesh handling.
Definition: mesh.hpp:592
snake
Definition: snake.hpp:83
mesh.hpp
Provides all the mesh tools used for the generation of 3D grids and geometries.
VertexNormal
std::tuple< coordvec, double > VertexNormal(const std::vector< double > &centre, const grid::coordlist &vecPts)
Calculates the vertex normal weighted by surface angle partitions.
Definition: mesh.cpp:747
meshprocessing.hpp
Tools for the mathematical processing of meshes.
postprocessing.hpp
Provide tecplot file formating for mesh and snake outputs.
triangulate.hpp
Provides a triangulation for snake and meshes.
warning.hpp
Provides the error and warning system used by the RSVS3D project.
RSVS3D_ERROR_NOTHROW
#define RSVS3D_ERROR_NOTHROW(M)
Generic rsvs warning.
Definition: warning.hpp:120
RSVS3D_ERROR_ARGUMENT
#define RSVS3D_ERROR_ARGUMENT(M)
Throw a invalid_argument.
Definition: warning.hpp:148