1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkEnSightWriter.cxx
5 
6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7   All rights reserved.
8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10      This software is distributed WITHOUT ANY WARRANTY; without even
11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12      PURPOSE.  See the above copyright notice for more information.
13 
14 =========================================================================*/
15 
16 /*----------------------------------------------------------------------------
17  Copyright (c) Sandia Corporation
18  See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
19 ----------------------------------------------------------------------------*/
20 
21 /* TODO
22  *
23  *
24  * check that data was written
25  *
26  */
27 
28 #include "vtkEnSightWriter.h"
29 
30 #include "vtkToolkits.h" // for VTK_USE_PARALLEL
31 #ifdef VTK_USE_PARALLEL
32 # include "vtkMultiProcessController.h"
33 #endif
34 
35 #include "vtkBitArray.h"
36 #include "vtkByteSwap.h"
37 #include "vtkCellArray.h"
38 #include "vtkCellData.h"
39 #include "vtkCharArray.h"
40 #include "vtkDataSet.h"
41 #include "vtkDoubleArray.h"
42 #include "vtkErrorCode.h"
43 #include "vtkFieldData.h"
44 #include "vtkFloatArray.h"
45 #include "vtkInformation.h"
46 #include "vtkIntArray.h"
47 #include "vtkLongArray.h"
48 #include "vtkLookupTable.h"
49 #include "vtkMath.h"
50 #include "vtkObjectFactory.h"
51 #include "vtkPointData.h"
52 #include "vtkPoints.h"
53 #include "vtkShortArray.h"
54 #include "vtkStreamingDemandDrivenPipeline.h"
55 #include "vtkUnsignedCharArray.h"
56 #include "vtkUnsignedIntArray.h"
57 #include "vtkUnsignedLongArray.h"
58 #include "vtkUnsignedShortArray.h"
59 
60 #include "vtkPriorityQueue.h"
61 
62 #include "vtkUnstructuredGrid.h"
63 
64 #include <cerrno>
65 #include <cmath>
66 #include <cctype>
67 
68 #include <vector>
69 #include <list>
70 #include <map>
71 #include <algorithm>
72 
73 // If we are building against a slightly older VTK version,
74 // these cell types are not defined, and won't occur in the input
75 
76 #ifndef VTK_QUADRATIC_WEDGE
77 # define VTK_QUADRATIC_WEDGE      26
78 # define VTK_QUADRATIC_PYRAMID    27
79 #endif
80 
81 // this undef is required on the hp. vtkMutexLock ends up including
82 // /usr/include/dce/cma_ux.h which has the gall to #define write as cma_write
83 
84 #ifdef write
85 # undef write
86 #endif
87 
88 //----------------------------------------------------------------------------
89 vtkStandardNewMacro(vtkEnSightWriter);
90 
91 //----------------------------------------------------------------------------
92 // Created object with no filename and timestep 0
vtkEnSightWriter()93 vtkEnSightWriter::vtkEnSightWriter()
94 {
95 
96   this->BaseName = nullptr;
97   this->FileName = nullptr;
98   this->TimeStep = 0;
99   this->Path=nullptr;
100   this->GhostLevelMultiplier=10000;
101   this->GhostLevel = 0;
102   this->TransientGeometry=false;
103   this->ProcessNumber=0;
104   this->NumberOfProcesses=1;
105   this->NumberOfBlocks=0;
106   this->BlockIDs=nullptr;
107   this->TmpInput = nullptr;
108 }
109 
110 //----------------------------------------------------------------------------
~vtkEnSightWriter()111 vtkEnSightWriter::~vtkEnSightWriter()
112 {
113   this->SetBaseName(nullptr);
114   this->SetFileName(nullptr);
115   this->SetPath(nullptr);
116 }
117 
118 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)119 void vtkEnSightWriter::PrintSelf(ostream& os, vtkIndent indent)
120 {
121   this->Superclass::PrintSelf(os,indent);
122 
123   os << indent << "FileName: "
124     << (this->FileName ? this->FileName : "(none)") << "\n";
125   os << indent << "Path: "
126     << (this->Path ? this->Path : "(none)") << "\n";
127   os << indent << "BaseName: "
128     << (this->BaseName ? this->BaseName : "(none)") << "\n";
129 
130   os << indent << "TimeStep: " << this->TimeStep << "\n";
131   os << indent << "TransientGeometry: " << this->TransientGeometry << "\n";
132   os << indent << "ProcessNumber: " << this->ProcessNumber << endl;
133   os << indent << "NumberOfProcesses: " << this->NumberOfProcesses << endl;
134   os << indent << "NumberOfBlocks: " << this->NumberOfBlocks << endl;
135   os << indent << "BlockIDs: " << this->BlockIDs << endl;
136   os << indent << "GhostLevel: " << this->GhostLevel << endl;
137 }
138 
FillInputPortInformation(int vtkNotUsed (port),vtkInformation * info)139 int vtkEnSightWriter::FillInputPortInformation( int vtkNotUsed(port), vtkInformation* info )
140 {
141   info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
142   return 1;
143 }
144 
145 //----------------------------------------------------------------------------
146 // Specify the input data or filter.
SetInputData(vtkUnstructuredGrid * input)147 void vtkEnSightWriter::SetInputData(vtkUnstructuredGrid *input)
148 {
149   this->SetInputDataInternal(0, input);
150 }
151 
152 
153 //----------------------------------------------------------------------------
154 // Specify the input data or filter.
GetInput()155 vtkUnstructuredGrid *vtkEnSightWriter::GetInput()
156 {
157   if (this->GetNumberOfInputConnections(0) < 1)
158   {
159     return nullptr;
160   }
161   else if (this->TmpInput)
162   {
163     return this->TmpInput;
164   }
165   else
166   {
167     return (vtkUnstructuredGrid *)(this->Superclass::GetInput());
168   }
169 }
170 
171 //----------------------------------------------------------------------------
WriteData()172 void vtkEnSightWriter::WriteData()
173 {
174   int i;
175   unsigned int ui;
176   int blockCount=0;
177   std::list<int>::iterator iter;
178 
179   if (this->TmpInput)
180   {
181     this->TmpInput->Delete();
182     this->TmpInput = nullptr;
183   }
184 
185   //figure out process ID
186 
187   this->ProcessNumber = 0;
188   this->NumberOfProcesses = 1;
189 
190 #ifdef VTK_USE_PARALLEL
191   vtkMultiProcessController *c = vtkMultiProcessController::GetGlobalController();
192 
193   if (c != nullptr)
194   {
195     this->ProcessNumber=c->GetLocalProcessId();
196     this->NumberOfProcesses = c->GetNumberOfProcesses();
197   }
198 #endif
199 
200   vtkUnstructuredGrid *input=this->GetInput();
201   vtkInformation* inInfo = this->GetInputInformation();
202 
203   if (this->GhostLevel >
204       inInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()))
205   {
206     // re-execute pipeline if necessary to obtain ghost cells
207 
208     this->GetInputAlgorithm()->UpdateInformation();
209     inInfo->Set(
210       vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
211       this->GhostLevel);
212     this->GetInputAlgorithm()->Update();
213   }
214 
215   //get the BlockID Cell Array
216   vtkDataArray *BlockData=input->GetCellData()->GetScalars("BlockId");
217 
218   if (BlockData==nullptr || strcmp(BlockData->GetName(),"BlockId"))
219   {
220     BlockData=nullptr;
221   }
222 
223   this->ComputeNames();
224 
225   if (!this->BaseName)
226   {
227     vtkErrorMacro("A FileName or Path/BaseName must be specified.");
228     return;
229   }
230 
231   this->SanitizeFileName(this->BaseName);
232 
233   char** blockNames=nullptr;
234   int * elementIDs=nullptr;
235   char charBuffer[1024];
236   char fileBuffer[512];
237   snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.%05d.geo",
238     this->Path,this->BaseName,this->ProcessNumber,
239     this->TimeStep);
240 
241   //open the geometry file
242   //only if timestep 0 and not transient geometry or transient geometry
243   FILE *fd=nullptr;
244   if (this->ShouldWriteGeometry())
245   {
246     if (!(fd=OpenFile(charBuffer)))
247       return;
248   }
249 
250   //Get the FILE's for Point Data Fields
251   std::vector<FILE*> pointArrayFiles;
252   int NumPointArrays=input->GetPointData()->GetNumberOfArrays();
253   for (i=0;i<NumPointArrays;i++)
254   {
255     strcpy(fileBuffer,input->GetPointData()->GetArray(i)->GetName());
256     this->SanitizeFileName(fileBuffer);
257     snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.%05d_n.%s",this->Path,this->BaseName,this->ProcessNumber,
258       this->TimeStep,fileBuffer);
259     FILE* ftemp=OpenFile(charBuffer);
260     if (!ftemp)
261     {
262       fclose(fd);
263       return;
264     }
265     pointArrayFiles.push_back(ftemp);
266 
267     //write the description line to the file
268     WriteStringToFile(fileBuffer,ftemp);
269   }
270 
271   //Get the FILE's for Cell Data Fields
272   std::vector<FILE*> cellArrayFiles;
273   int NumCellArrays=input->GetCellData()->GetNumberOfArrays();
274   for (i=0;i<NumCellArrays;i++)
275   {
276 
277     strcpy(fileBuffer,input->GetCellData()->GetArray(i)->GetName());
278     this->SanitizeFileName(fileBuffer);
279     snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.%05d_c.%s",this->Path,this->BaseName,this->ProcessNumber,
280       this->TimeStep,fileBuffer);
281     FILE* ftemp=OpenFile(charBuffer);
282     if (!ftemp)
283     {
284       fclose(fd);
285       return;
286     }
287     cellArrayFiles.push_back(ftemp);
288 
289     //write the description line to the file
290     WriteStringToFile(fileBuffer,ftemp);
291   }
292 
293   //write the header information
294   if (this->ShouldWriteGeometry())
295   {
296     this->WriteStringToFile("C Binary",fd);
297     this->WriteStringToFile("Written by VTK EnSight Writer",fd);
298     //if (this->Title)
299     // this->WriteStringToFile(this->Title,fd);
300     //else
301     this->WriteStringToFile("No Title was Specified",fd);
302     //we will specify node and element ID's
303     this->WriteStringToFile("node id given\n",fd);
304     this->WriteStringToFile("element id given\n",fd);
305   }
306 
307   //get the Ghost Cell Array if it exists
308   vtkDataArray *GhostData=input->GetCellData()->GetScalars(vtkDataSetAttributes::GhostArrayName());
309   //if the strings are not the same then we did not get the ghostData array
310   if (GhostData==nullptr || strcmp(GhostData->GetName(), vtkDataSetAttributes::GhostArrayName()))
311   {
312     GhostData=nullptr;
313   }
314 
315 
316   //data structure to get all the cells for a certain part
317   //basically sort by part# and cell type
318   std::map <int, std::vector <int> > CellsByPart;
319 
320   //just a list of part numbers
321   std::list<int> partNumbers;
322 
323   //get all the part numbers in the unstructured grid and sort the cells
324   //by part number
325   for (i=0;i<input->GetNumberOfCells();i++)
326   {
327     int key=1;
328     if (BlockData)
329       key=(int)(BlockData->GetTuple(i)[0]);
330     else
331       cout << "No BlockID was found\n";
332     if (CellsByPart.count(key)==0)
333     {
334       CellsByPart[key]=std::vector < int >() ;
335     }
336     CellsByPart[key].push_back(i);
337     partNumbers.push_back(key);
338 
339   }
340 
341   //remove the duplicates from the partNumbers
342   partNumbers.sort();
343   partNumbers.unique();
344 
345   //write out each part
346   for (iter=partNumbers.begin();iter!=partNumbers.end();++iter)
347   {
348     unsigned int j;
349     std::list<int>::iterator iter2;
350     int part=*iter;
351 
352     //write the part Header
353     if (this->ShouldWriteGeometry())
354     {
355       blockCount+=1;
356       this->WriteStringToFile("part",fd);
357       this->WriteIntToFile(part,fd);
358       //cout << "part is " << part << endl;
359       this->WriteStringToFile("VTK Part",fd);
360       this->WriteStringToFile("coordinates",fd);
361     }
362 
363     //write the part header for data files
364     for (j=0;j<pointArrayFiles.size();j++)
365     {
366       this->WriteStringToFile("part",pointArrayFiles[j]);
367       this->WriteIntToFile(part,pointArrayFiles[j]);
368       this->WriteStringToFile("coordinates",pointArrayFiles[j]);
369     }
370     for (j=0;j<cellArrayFiles.size();j++)
371     {
372       this->WriteStringToFile("part",cellArrayFiles[j]);
373       this->WriteIntToFile(part,cellArrayFiles[j]);
374     }
375 
376     //list of VTK Node Indices per part
377     std::list<int> NodesPerPart;
378 
379     //map that goes from NodeID to the order, used for element connectivity
380     std::map<int, int> NodeIdToOrder;
381 
382     //get a list of all the nodes used for a particular part
383     for (j=0;j<CellsByPart[part].size();j++)
384     {
385       vtkCell* cell=input->GetCell(CellsByPart[part][j]);
386       vtkIdList* points=cell->GetPointIds();
387       for (int k=0;k<points->GetNumberOfIds();k++)
388       {
389         NodesPerPart.push_back(points->GetId(k));
390       }
391     }
392 
393     //remove the duplicate Node ID's
394     NodesPerPart.sort();
395     NodesPerPart.unique();
396 
397     //write the number of nodes
398     if (this->ShouldWriteGeometry())
399     {
400       this->WriteIntToFile(static_cast<int>(NodesPerPart.size()),fd);
401 
402 
403       //write the Node ID's to the file
404       //also set up the NodeID->order map
405       int NodeCount=0;
406       for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
407       {
408         this->WriteIntToFile(*iter2,fd);
409         NodeIdToOrder[*iter2]=NodeCount+1;
410         NodeCount++;
411       }
412 
413 
414       //EnSight format requires all the X's, then all the Y's, then all the Z's
415       //write the X Coordinates
416 
417       vtkPoints* inputPoints=input->GetPoints();
418       for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
419       {
420         this->WriteFloatToFile((float)(inputPoints->GetPoint(*iter2)[0]),fd);
421       }
422       for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
423       {
424         this->WriteFloatToFile((float)(inputPoints->GetPoint(*iter2)[1]),fd);
425       }
426       for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
427       {
428         this->WriteFloatToFile((float)(inputPoints->GetPoint(*iter2)[2]),fd);
429       }
430 
431     }
432 
433     //write the Node Data for this part
434     for (j=0;j<pointArrayFiles.size();j++)
435     {
436       vtkDataArray* DataArray=input->GetPointData()->GetArray(j);
437       //figure out what type of data it is
438       int DataSize=DataArray->GetNumberOfComponents();
439 
440       for (int CurrentDimension=0;
441         CurrentDimension<DataSize;
442         CurrentDimension++)
443       {
444         for (std::list<int>::iterator k=NodesPerPart.begin();
445           k!=NodesPerPart.end();++k)
446         {
447           this->WriteFloatToFile((float)
448             (DataArray->
449              GetTuple(*k)[CurrentDimension]),
450             pointArrayFiles[j]);
451         }
452       }
453     }
454 
455 
456 
457     //now we need to sort the cell list by element type
458     //map is indexed by cell type has a vector of cell ID's
459     std::map<int, std::vector<int> > CellsByElement;
460     for (j=0;j<CellsByPart[part].size();j++)
461     {
462       int CellType=input->GetCell(CellsByPart[part][j])->GetCellType();
463       int ghostLevel=0;
464       if (GhostData)
465       {
466         ghostLevel=(int)(GhostData->GetTuple(CellsByPart[part][j])[0]);
467         if (ghostLevel & vtkDataSetAttributes::DUPLICATECELL)
468         {
469           ghostLevel=1;
470         }
471       }
472       //we want to sort out the ghost cells from the normal cells
473       //so the element type will be ghostMultiplier*ghostLevel+elementType
474       CellType+=ghostLevel*GhostLevelMultiplier;
475       if (CellsByElement.count(CellType)==0)
476       {
477         CellsByElement[CellType]=std::vector < int >() ;
478       }
479       CellsByElement[CellType].push_back(CellsByPart[part][j]);
480     }
481 
482     //now we need to go through each element type that EnSight understands
483     std::vector<int> elementTypes;
484 
485     //list the types that EnSight understands
486     //the noticeable absences are the ones without a fixed number of Nodes
487     //for the ghost cell types
488     elementTypes.push_back(VTK_VERTEX);
489     elementTypes.push_back(VTK_LINE);
490     elementTypes.push_back(VTK_TRIANGLE);
491     elementTypes.push_back(VTK_QUAD);
492     elementTypes.push_back(VTK_POLYGON);
493     elementTypes.push_back(VTK_TETRA);
494     elementTypes.push_back(VTK_HEXAHEDRON);
495     elementTypes.push_back(VTK_WEDGE);
496     elementTypes.push_back(VTK_PYRAMID);
497     elementTypes.push_back(VTK_CONVEX_POINT_SET);
498     elementTypes.push_back(VTK_QUADRATIC_EDGE);
499     elementTypes.push_back(VTK_QUADRATIC_TRIANGLE);
500     elementTypes.push_back(VTK_QUADRATIC_QUAD);
501     elementTypes.push_back(VTK_QUADRATIC_TETRA);
502     elementTypes.push_back(VTK_QUADRATIC_HEXAHEDRON);
503     elementTypes.push_back(VTK_QUADRATIC_WEDGE);
504     elementTypes.push_back(VTK_QUADRATIC_PYRAMID);
505     elementTypes.push_back(this->GhostLevelMultiplier+VTK_VERTEX);
506     elementTypes.push_back(this->GhostLevelMultiplier+VTK_LINE);
507     elementTypes.push_back(this->GhostLevelMultiplier+VTK_TRIANGLE);
508     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUAD);
509     elementTypes.push_back(this->GhostLevelMultiplier+VTK_POLYGON);
510     elementTypes.push_back(this->GhostLevelMultiplier+VTK_TETRA);
511     elementTypes.push_back(this->GhostLevelMultiplier+VTK_HEXAHEDRON);
512     elementTypes.push_back(this->GhostLevelMultiplier+VTK_WEDGE);
513     elementTypes.push_back(this->GhostLevelMultiplier+VTK_PYRAMID);
514     elementTypes.push_back(this->GhostLevelMultiplier+VTK_CONVEX_POINT_SET);
515     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_EDGE);
516     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_TRIANGLE);
517     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_QUAD);
518     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_TETRA);
519     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_HEXAHEDRON);
520     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_WEDGE);
521     elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_PYRAMID);
522 
523     //write out each type of element
524     if (this->ShouldWriteGeometry())
525     {
526       for (j=0;j<elementTypes.size();j++)
527       {
528         unsigned int k;
529         int elementType=elementTypes[j];
530         if (CellsByElement.count(elementType)>0)
531         {
532           //switch on element type to write correct type to file
533           this->WriteElementTypeToFile(elementType,fd);
534 
535           //number of elements
536           this->WriteIntToFile(
537             static_cast<int>(CellsByElement[elementType].size()),fd);
538 
539           //element ID's
540           for (k=0;k<CellsByElement[elementType].size();k++)
541           {
542             int CellId=CellsByElement[elementType][k];
543             this->WriteIntToFile(CellId,fd);
544           }
545 
546           //element conenctivity information
547           for (k=0;k<CellsByElement[elementType].size();k++)
548           {
549             int CellId=CellsByElement[elementType][k];
550             vtkIdList *PointIds=input->GetCell(CellId)->GetPointIds();
551             for (int m=0;m<PointIds->GetNumberOfIds();m++)
552             {
553               int PointId=PointIds->GetId(m);
554               this->WriteIntToFile(NodeIdToOrder[PointId],fd);
555             }
556           }
557         }
558       }
559     }
560 
561     //write the Cell Data for this part
562     for (j=0;j<cellArrayFiles.size();j++)
563     {
564       vtkDataArray* DataArray=input->GetCellData()->GetArray(j);
565       //figure out what type of data it is
566       int DataSize=DataArray->GetNumberOfComponents();
567 
568       for (unsigned int k=0;k<elementTypes.size();k++)
569       {
570         if (!CellsByElement[elementTypes[k]].empty())
571         {
572           this->WriteElementTypeToFile(elementTypes[k],
573             cellArrayFiles[j]);
574           for (unsigned int m=0;m<CellsByElement[elementTypes[k]].size();m++)
575           {
576             for ( int CurrentDimension = 0; CurrentDimension < DataSize;
577                      CurrentDimension ++ )
578             {
579               this->WriteFloatToFile
580                 (  (float)
581                    (   DataArray
582                        ->GetTuple(  CellsByElement[ elementTypes[k] ][m]  )
583                          [CurrentDimension]
584                    ),
585                    cellArrayFiles[j]
586                 );
587             }
588           }
589         }
590       }
591     }
592   }
593 
594   //now write the empty blocks
595   //use the block list in the exodus model if it exists, otherwise
596   //use the BlockID list if that exists.
597 
598   if (this->BlockIDs)
599   {
600     elementIDs=this->BlockIDs;
601   }
602 
603   if (elementIDs)
604   {
605     //cout << "have " << this->NumberOfBlocks << " blocks " << endl;
606     for (i=0;i<this->NumberOfBlocks;i++)
607     {
608       unsigned int j;
609       //figure out if the part was already written
610       int part=elementIDs[i];
611       if (   std::find(partNumbers.begin(), partNumbers.end(), part)
612         == partNumbers.end() )
613       {
614         //no information about the part was written to the output files
615         //so write some empty information
616         if (this->ShouldWriteGeometry())
617         {
618           blockCount+=1;
619           this->WriteStringToFile("part",fd);
620           this->WriteIntToFile(part,fd);
621 
622           int exodusIndex =
623             this->GetExodusModelIndex(elementIDs,this->NumberOfBlocks,part);
624 
625           if (exodusIndex!=-1 && blockNames)
626           {
627             snprintf(charBuffer,sizeof(charBuffer),"Exodus-%s-%d",blockNames[exodusIndex],part);
628             this->WriteStringToFile(charBuffer,fd);
629           }
630           else
631           {
632             this->WriteStringToFile("VTK Part",fd);
633           }
634         }
635 
636         //write the part header for data files
637         for (j=0;j<pointArrayFiles.size();j++)
638         {
639           this->WriteStringToFile("part",pointArrayFiles[j]);
640           this->WriteIntToFile(part,pointArrayFiles[j]);
641         }
642         for (j=0;j<cellArrayFiles.size();j++)
643         {
644           this->WriteStringToFile("part",cellArrayFiles[j]);
645           this->WriteIntToFile(part,cellArrayFiles[j]);
646         }
647       }
648     }
649   }
650   //cout << "wrote " << blockCount << "parts\n";
651   if (this->TmpInput)
652   {
653     this->TmpInput->Delete();
654     this->TmpInput = nullptr;
655   }
656 
657   //close all the files
658   if (fd)
659   {
660     fclose(fd);
661   }
662 
663   for (ui=0;ui<cellArrayFiles.size();ui++)
664   {
665     fclose(cellArrayFiles[ui]);
666   }
667   for (ui=0;ui<pointArrayFiles.size();ui++)
668   {
669     fclose(pointArrayFiles[ui]);
670   }
671 
672 }
673 
674 //----------------------------------------------------------------------------
WriteCaseFile(int TotalTimeSteps)675 void vtkEnSightWriter::WriteCaseFile(int TotalTimeSteps)
676 {
677 
678   vtkUnstructuredGrid *input=this->GetInput();
679   int i;
680 
681   this->ComputeNames();
682 
683   if (!this->BaseName)
684   {
685     vtkErrorMacro("A FileName or Path/BaseName must be specified.");
686     return;
687   }
688 
689   char charBuffer[1024];
690   snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.case",this->Path,this->BaseName,this->ProcessNumber);
691 
692   //open the geometry file
693   FILE *fd=nullptr;
694   if (!(fd=OpenFile(charBuffer)))
695   {
696     return;
697   }
698 
699   this->WriteTerminatedStringToFile("FORMAT\n",fd);
700   this->WriteTerminatedStringToFile("type: ensight gold\n\n",fd);
701   this->WriteTerminatedStringToFile("\nGEOMETRY\n",fd);
702 
703   //write the geometry file
704   if (!this->TransientGeometry)
705   {
706     snprintf(charBuffer,sizeof(charBuffer),"model: %s.%d.00000.geo\n",this->BaseName,this->ProcessNumber);
707     this->WriteTerminatedStringToFile(charBuffer,fd);
708   }
709   else
710   {
711     snprintf(charBuffer,sizeof(charBuffer),"model: 1 %s.%d.*****.geo\n",this->BaseName,this->ProcessNumber);
712     this->WriteTerminatedStringToFile(charBuffer,fd);
713   }
714 
715   this->WriteTerminatedStringToFile("\nVARIABLE\n",fd);
716 
717   char fileBuffer[256];
718 
719   //write the Node variable files
720   for (i=0;i<input->GetPointData()->GetNumberOfArrays();i++)
721   {
722 
723     strcpy(fileBuffer,input->GetPointData()->GetArray(i)->GetName());
724     // skip arrays that were not written
725     if (strcmp(fileBuffer,"GlobalElementId")==0)
726     {
727       continue;
728     }
729     if (strcmp(fileBuffer,"GlobalNodeId")==0)
730     {
731       continue;
732     }
733     if (strcmp(fileBuffer,"BlockId")==0)
734     {
735       continue;
736     }
737     this->SanitizeFileName(fileBuffer);
738     //figure out what kind of data it is
739     char SmallBuffer[16];
740     switch(input->GetPointData()->GetArray(i)->GetNumberOfComponents())
741     {
742     case(1):
743       strcpy(SmallBuffer,"scalar");
744       break;
745     case(3):
746       strcpy(SmallBuffer,"vector");
747       break;
748     case(6):
749       strcpy(SmallBuffer,"tensor");
750       break;
751     case(9):
752       strcpy(SmallBuffer,"tensor9");
753       break;
754     }
755     if (TotalTimeSteps<=1)
756     {
757       snprintf(charBuffer,sizeof(charBuffer),"%s per node: %s_n %s.%d.00000_n.%s\n",
758         SmallBuffer,
759         fileBuffer,
760         this->BaseName,
761         this->ProcessNumber,
762         fileBuffer);
763     }
764     else
765     {
766       snprintf(charBuffer,sizeof(charBuffer),"%s per node: 1 %s_n %s.%d.*****_n.%s\n",
767         SmallBuffer,
768         fileBuffer,
769         this->BaseName,
770         this->ProcessNumber,
771         fileBuffer);
772 
773 
774     }
775     this->WriteTerminatedStringToFile(charBuffer,fd);
776   }
777 
778   //write the cell variable files
779   for (i=0;i<input->GetCellData()->GetNumberOfArrays();i++)
780   {
781     //figure out what kind of data it is
782     char SmallBuffer[16];
783 
784     strcpy(fileBuffer,input->GetCellData()->GetArray(i)->GetName());
785     // skip arrays that were not written
786     if (strcmp(fileBuffer,"GlobalElementId")==0)
787     {
788       continue;
789     }
790     if (strcmp(fileBuffer,"GlobalNodeId")==0)
791     {
792       continue;
793     }
794     if (strcmp(fileBuffer,"BlockId")==0)
795     {
796       continue;
797     }
798     this->SanitizeFileName(fileBuffer);
799     switch(input->GetCellData()->GetArray(i)->GetNumberOfComponents())
800     {
801     case(1):
802       strcpy(SmallBuffer,"scalar");
803       break;
804     case(3):
805       strcpy(SmallBuffer,"vector");
806       break;
807     case(6):
808       strcpy(SmallBuffer,"tensor");
809       break;
810     case(9):
811       strcpy(SmallBuffer,"tensor9");
812       break;
813     }
814     if (TotalTimeSteps<=1)
815     {
816       snprintf(charBuffer,sizeof(charBuffer),"%s per element: %s_c %s.%d.00000_c.%s\n",
817         SmallBuffer,
818         fileBuffer,
819         this->BaseName,
820         this->ProcessNumber,
821         fileBuffer);
822     }
823     else
824     {
825       snprintf(charBuffer,sizeof(charBuffer),"%s per element: 1 %s_c %s.%d.*****_c.%s\n",
826         SmallBuffer,
827         fileBuffer,
828         this->BaseName,
829         this->ProcessNumber,
830         fileBuffer);
831     }
832     this->WriteTerminatedStringToFile(charBuffer,fd);
833   }
834 
835   //write time information if we have multiple timesteps
836   if (TotalTimeSteps>1)
837   {
838     this->WriteTerminatedStringToFile("\nTIME\n",fd);
839     this->WriteTerminatedStringToFile("time set: 1\n",fd);
840     snprintf(charBuffer,sizeof(charBuffer),"number of steps: %d\n",TotalTimeSteps);
841     this->WriteTerminatedStringToFile(charBuffer,fd);
842     this->WriteTerminatedStringToFile("filename start number: 00000\n",fd);
843     this->WriteTerminatedStringToFile("filename increment: 00001\n",fd);
844     this->WriteTerminatedStringToFile("time values: \n",fd);
845     for (i=0;i<TotalTimeSteps;i++)
846     {
847       double timestep=i;
848 
849       snprintf(charBuffer,sizeof(charBuffer),"%f ",timestep);
850       this->WriteTerminatedStringToFile(charBuffer,fd);
851       if (i%6==0 && i>0)
852       {
853         this->WriteTerminatedStringToFile("\n",fd);
854       }
855     }
856   }
857 
858 
859 }
860 
861 //----------------------------------------------------------------------------
WriteSOSCaseFile(int numProcs)862 void vtkEnSightWriter::WriteSOSCaseFile(int numProcs)
863 {
864   this->ComputeNames();
865 
866   if (!this->BaseName)
867   {
868     vtkErrorMacro("A FileName or Path/BaseName must be specified.");
869     return;
870   }
871 
872   this->SanitizeFileName(this->BaseName);
873 
874   char charBuffer[512];
875   snprintf(charBuffer,sizeof(charBuffer),"%s/%s.case.sos",this->Path,this->BaseName);
876 
877   FILE *fd=nullptr;
878   if (!(fd=OpenFile(charBuffer)))
879     return;
880 
881   this->WriteTerminatedStringToFile("FORMAT\n",fd);
882   this->WriteTerminatedStringToFile("type: master_server gold\n\n",fd);
883 
884   this->WriteTerminatedStringToFile("SERVERS\n",fd);
885   snprintf(charBuffer,sizeof(charBuffer),"number of servers: %d\n\n",numProcs);
886   this->WriteTerminatedStringToFile(charBuffer,fd);
887 
888   //write the servers section with placeholders for the ensight server
889   //location and server name
890   int i=0;
891   for (i=0;i<numProcs;i++)
892   {
893     snprintf(charBuffer,sizeof(charBuffer), "#Server %d\n",i);
894     this->WriteTerminatedStringToFile(charBuffer,fd);
895     this->WriteTerminatedStringToFile("#-------\n",fd);
896     snprintf(charBuffer,sizeof(charBuffer), "machine id: MID%05d\n",i);
897     this->WriteTerminatedStringToFile(charBuffer,fd);
898 
899     this->WriteTerminatedStringToFile("executable: MEX\n",fd);
900     snprintf(charBuffer,sizeof(charBuffer), "data_path: %s\n",this->Path);
901     this->WriteTerminatedStringToFile(charBuffer,fd);
902     snprintf(charBuffer,sizeof(charBuffer),"casefile: %s.%d.case\n\n",this->BaseName,i);
903     this->WriteTerminatedStringToFile(charBuffer,fd);
904   }
905 
906 
907 }
908 
909 //----------------------------------------------------------------------------
WriteStringToFile(const char * cstring,FILE * file)910 void vtkEnSightWriter::WriteStringToFile(const char* cstring, FILE* file)
911 {
912   char cbuffer[81];
913   unsigned long cstringLength = static_cast<unsigned long>(strlen(cstring));
914   memcpy(cbuffer,cstring,vtkMath::Min(cstringLength,80ul));
915   for (int i = cstringLength; i < 81; ++i)
916   {
917     cbuffer[i] = '\0';
918   }
919 
920   // Write a constant 80 bytes to the file.
921   fwrite(cbuffer, sizeof(char),80,file);
922 }
923 
924 //----------------------------------------------------------------------------
WriteTerminatedStringToFile(const char * cstring,FILE * file)925 void vtkEnSightWriter::WriteTerminatedStringToFile(const char* cstring, FILE* file)
926 {
927   fwrite(cstring, sizeof(char),std::min(strlen(cstring), static_cast<size_t>(512)),file);
928 }
929 
930 //----------------------------------------------------------------------------
WriteIntToFile(const int i,FILE * file)931 void vtkEnSightWriter::WriteIntToFile(const int i,FILE* file)
932 {
933   fwrite(&i, sizeof(int),1,file);
934 }
935 
936 //----------------------------------------------------------------------------
WriteFloatToFile(const float f,FILE * file)937 void vtkEnSightWriter::WriteFloatToFile(const float f,FILE* file)
938 {
939   fwrite(&f, sizeof(float),1,file);
940 }
941 
942 //----------------------------------------------------------------------------
WriteElementTypeToFile(int elementType,FILE * fd)943 void vtkEnSightWriter::WriteElementTypeToFile(int elementType,FILE* fd)
944 {
945   int ghostLevel=elementType/GhostLevelMultiplier;
946   elementType=elementType%GhostLevelMultiplier;
947   if (ghostLevel==0)
948   {
949     switch(elementType)
950     {
951     case(VTK_VERTEX):
952       this->WriteStringToFile("point",fd);
953       break;
954     case(VTK_LINE):
955       this->WriteStringToFile("bar2",fd);
956       break;
957     case(VTK_TRIANGLE):
958       this->WriteStringToFile("tria3",fd);
959       break;
960     case(VTK_QUAD):
961       this->WriteStringToFile("quad4",fd);
962       break;
963     case(VTK_POLYGON):
964       this->WriteStringToFile("nsided",fd);
965       break;
966     case(VTK_TETRA):
967       this->WriteStringToFile("tetra4",fd);
968       break;
969     case(VTK_HEXAHEDRON):
970       this->WriteStringToFile("hexa8",fd);
971       break;
972     case(VTK_WEDGE):
973       this->WriteStringToFile("penta6",fd);
974       break;
975     case(VTK_PYRAMID):
976       this->WriteStringToFile("pyramid5",fd);
977       break;
978     case(VTK_CONVEX_POINT_SET):
979       this->WriteStringToFile("nfaced",fd);
980       break;
981     case(VTK_QUADRATIC_EDGE):
982       this->WriteStringToFile("bar3",fd);
983       break;
984     case(VTK_QUADRATIC_TRIANGLE):
985       this->WriteStringToFile("tria6",fd);
986       break;
987     case(VTK_QUADRATIC_QUAD):
988       this->WriteStringToFile("quad8",fd);
989       break;
990     case(VTK_QUADRATIC_TETRA):
991       this->WriteStringToFile("tetra10",fd);
992       break;
993     case(VTK_QUADRATIC_HEXAHEDRON):
994       this->WriteStringToFile("hexa20",fd);
995       break;
996     case(VTK_QUADRATIC_WEDGE):
997       this->WriteStringToFile("penta15",fd);
998       break;
999     case(VTK_QUADRATIC_PYRAMID):
1000       this->WriteStringToFile("pyramid13",fd);
1001       break;
1002     }
1003   }
1004   else
1005   {
1006     switch(elementType)
1007     {
1008     case(VTK_VERTEX):
1009       this->WriteStringToFile("g_point",fd);
1010       break;
1011     case(VTK_LINE):
1012       this->WriteStringToFile("g_bar2",fd);
1013       break;
1014     case(VTK_TRIANGLE):
1015       this->WriteStringToFile("g_tria3",fd);
1016       break;
1017     case(VTK_QUAD):
1018       this->WriteStringToFile("g_quad4",fd);
1019       break;
1020     case(VTK_POLYGON):
1021       this->WriteStringToFile("g_nsided",fd);
1022       break;
1023     case(VTK_TETRA):
1024       this->WriteStringToFile("g_tetra4",fd);
1025       break;
1026     case(VTK_HEXAHEDRON):
1027       this->WriteStringToFile("g_hexa8",fd);
1028       break;
1029     case(VTK_WEDGE):
1030       this->WriteStringToFile("g_penta6",fd);
1031       break;
1032     case(VTK_PYRAMID):
1033       this->WriteStringToFile("g_pyramid5",fd);
1034       break;
1035     case(VTK_CONVEX_POINT_SET):
1036       this->WriteStringToFile("g_nfaced",fd);
1037       break;
1038     case(VTK_QUADRATIC_EDGE):
1039       this->WriteStringToFile("g_bar3",fd);
1040       break;
1041     case(VTK_QUADRATIC_TRIANGLE):
1042       this->WriteStringToFile("g_tria6",fd);
1043       break;
1044     case(VTK_QUADRATIC_QUAD):
1045       this->WriteStringToFile("g_quad8",fd);
1046       break;
1047     case(VTK_QUADRATIC_TETRA):
1048       this->WriteStringToFile("g_tetra10",fd);
1049       break;
1050     case(VTK_QUADRATIC_HEXAHEDRON):
1051       this->WriteStringToFile("g_hexa20",fd);
1052       break;
1053     case(VTK_QUADRATIC_WEDGE):
1054       this->WriteStringToFile("g_penta15",fd);
1055       break;
1056     case(VTK_QUADRATIC_PYRAMID):
1057       this->WriteStringToFile("g_pyramid13",fd);
1058       break;
1059     }
1060   }
1061 }
1062 
1063 //----------------------------------------------------------------------------
ShouldWriteGeometry()1064 bool vtkEnSightWriter::ShouldWriteGeometry()
1065 {
1066   return (this->TransientGeometry || (this->TimeStep==0));
1067 }
1068 
1069 //----------------------------------------------------------------------------
SanitizeFileName(char * name)1070 void vtkEnSightWriter::SanitizeFileName(char* name)
1071 {
1072 
1073   char buffer[512];
1074   unsigned int i;
1075   int BufferPosition=0;
1076   for (i=0;i<strlen(name);i++)
1077   {
1078     if (name[i]!='/')
1079     {
1080       buffer[BufferPosition]=name[i];
1081       BufferPosition++;
1082     }
1083   }
1084   buffer[BufferPosition]=0;
1085   for (i=0;i<strlen(buffer);i++)
1086   {
1087     name[i]=buffer[i];
1088   }
1089   name[strlen(buffer)]=0;
1090 
1091 }
1092 
1093 //----------------------------------------------------------------------------
OpenFile(char * name)1094 FILE* vtkEnSightWriter::OpenFile(char* name)
1095 {
1096   FILE * fd=fopen(name,"wb");
1097 
1098   if (fd == nullptr)
1099   {
1100     vtkErrorMacro("Error opening " << name
1101       << ": " << strerror(errno));
1102     return nullptr;
1103   }
1104   return fd;
1105 }
1106 
1107 //----------------------------------------------------------------------------
GetExodusModelIndex(int * elementArray,int numberElements,int partID)1108 int vtkEnSightWriter::GetExodusModelIndex(int *elementArray,int numberElements,int partID)
1109 {
1110   int i;
1111   for (i=0;i<numberElements;i++)
1112   {
1113     if (elementArray[i]==partID)
1114       return i;
1115   }
1116   return -1;
1117 }
1118 
1119 //----------------------------------------------------------------------------
DefaultNames()1120 void vtkEnSightWriter::DefaultNames()
1121 {
1122   char *path = new char [4];
1123   char *base = new char [20];
1124   strcpy(path, "./");
1125   strcpy(base, "EnSightWriter.out");
1126 
1127   this->SetPath(path);
1128   this->SetBaseName(base);
1129 }
1130 
1131 //----------------------------------------------------------------------------
ComputeNames()1132 void vtkEnSightWriter::ComputeNames()
1133 {
1134   if (this->Path && this->BaseName)
1135   {
1136     return;
1137   }
1138 
1139   if (!this->FileName)
1140   {
1141     this->DefaultNames();
1142     return;
1143   }
1144 
1145   // FileName = Path/BaseName.digits.digits
1146 
1147   char *path = nullptr;
1148   char *base = nullptr;
1149 
1150   char *f = this->FileName;
1151 
1152   while (!isgraph(*f)) f++;  // find first printable character
1153 
1154   if (!*f)
1155   {
1156     // FileName is garbage
1157     DefaultNames();
1158     return;
1159   }
1160 
1161   char *buf = new char [strlen(f) + 1];
1162   strcpy(buf, f);
1163 
1164   char *slash = strrchr(buf, '/');  // final slash
1165 
1166   if (slash)
1167   {
1168     *slash = 0;
1169     path = new char [strlen(buf) + 1];
1170     strcpy(path, buf);
1171     f = slash + 1;
1172   }
1173   else
1174   {
1175     path = new char [4];
1176     strcpy(path, "./");
1177 
1178     f = buf;
1179   }
1180 
1181   char *firstChar = f;
1182   while (*f && (*f != '.')) f++;
1183   *f = 0;
1184 
1185   base = new char [strlen(firstChar) + 1];
1186   strcpy(base, firstChar);
1187 
1188   this->SetPath(path);
1189   this->SetBaseName(base);
1190 
1191   delete [] buf;
1192 }
1193