1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: $RCSfile: vtkPhastaReader.cxx,v $
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 #include "vtkPhastaReader.h"
16
17 #include "vtkByteSwap.h"
18 #include "vtkCellType.h" //added for constants such as VTK_TETRA etc...
19 #include "vtkDataArray.h"
20 #include "vtkIntArray.h"
21 #include "vtkDoubleArray.h"
22 #include "vtkFloatArray.h"
23 #include "vtkInformation.h"
24 #include "vtkInformationVector.h"
25 #include "vtkObjectFactory.h"
26 #include "vtkPointData.h"
27 #include "vtkCellData.h"
28 #include "vtkPointSet.h"
29 #include "vtkSmartPointer.h"
30 #include "vtkUnstructuredGrid.h"
31 //change
32 #include "vtkStreamingDemandDrivenPipeline.h"
33 #include "vtkPVXMLElement.h"
34 #include "vtkPVXMLParser.h"
35
36 #include "vtkCellData.h"
37 #include "vtkFieldData.h"
38 #include "vtkMultiBlockDataSet.h"
39 #include "vtkMultiPieceDataSet.h"
40
41 //change end
42 vtkCxxRevisionMacro(vtkPhastaReader, "$Revision: 1.11 $");
43 vtkStandardNewMacro(vtkPhastaReader);
44
45 #include <vtkstd/map>
46 #include <vtkstd/vector>
47 #include <vtkstd/string>
48 #include <vtksys/ios/sstream>
49
50 //CHANGE////////////////////////////////////////////////////////////////
51 #include "rdtsc.h"
52 #define clockRate 2670000000.0
53 unsigned long long start, end;
54 //double opentime_total = 0.0;
55 int LAST_FILE_ID;
56
startTimer(unsigned long long * start)57 void startTimer(unsigned long long* start) {
58 *start = rdtsc();
59 }
60
endTimer(unsigned long long * end)61 void endTimer(unsigned long long* end) {
62 *end = rdtsc();
63 }
64
computeTime(unsigned long long * start,unsigned long long * end)65 void computeTime(unsigned long long* start, unsigned long long* end) {
66 double time = (double)((*end-*start)/clockRate);
67 opentime_total += time;
68 }
69
70
71 #define VERSION_INFO_HEADER_SIZE 8192
72 #define DB_HEADER_SIZE 1024
73 #define TWO_MEGABYTE 2097152
74 #define ENDIAN_TEST_NUMBER 12180 // Troy's Zip Code!!
75 #define MAX_PHASTA_FILES 64
76 #define MAX_PHASTA_FILE_NAME_LENGTH 1024
77 #define MAX_FIELDS_NUMBER 48
78 #define MAX_FIELDS_NAME_LENGTH 128
79 #define DefaultMHSize (4*1024*1024)
80 int MasterHeaderSize = DefaultMHSize;
81 int diff_endian = 0;
82 long long counter = 0;
83
84 enum PhastaIO_Errors
85 {
86 MAX_PHASTA_FILES_EXCEEDED = -1,
87 UNABLE_TO_OPEN_FILE = -2,
88 NOT_A_MPI_FILE = -3,
89 GPID_EXCEEDED = -4,
90 DATA_TYPE_ILLEGAL = -5,
91 };
92
93 //CHANGE END////////////////////////////////////////////////////////////
94 struct vtkPhastaReaderInternal
95 {
96 struct FieldInfo
97 {
98 int StartIndexInPhastaArray;
99 int NumberOfComponents;
100 int DataDependency; // 0-nodal, 1-elemental
101 vtkstd::string DataType; // "int" or "double"
102 vtkstd::string PhastaFieldTag;
103
FieldInfovtkPhastaReaderInternal::FieldInfo104 FieldInfo() : StartIndexInPhastaArray(-1), NumberOfComponents(-1), DataDependency(-1), DataType(""), PhastaFieldTag("")
105 {
106 }
107 };
108
109 typedef vtkstd::map<vtkstd::string, FieldInfo> FieldInfoMapType;
110 FieldInfoMapType FieldInfoMap;
111 };
112
113
114 // Begin of copy from phastaIO
115
116 //CHANGE////////////////////////////////////////////////////////////////
117
118 /***********************************************************************/
119 /***************** NEW PHASTA IO CODE STARTS HERE **********************/
120 /***********************************************************************/
121
122 int partID_counter;
123
124 typedef struct
125 {
126 bool Wrong_Endian; /* default to false */
127
128 char filename[MAX_PHASTA_FILE_NAME_LENGTH]; /* defafults to 1024 */
129 int nppp;
130 int nPPF;
131 int nFiles;
132 int nFields;
133 unsigned long long my_offset;
134
135 char * master_header;
136 double * double_chunk;
137 int * int_chunk;
138 double * read_double_chunk;
139 int * read_int_chunk;
140 unsigned long long **my_offset_table;
141 unsigned long long **my_read_table;
142 int field_count;
143 int part_count;
144 int read_field_count;
145 int read_part_count;
146 int GPid;
147 int start_id;
148 unsigned long long next_start_address;
149
150 int myrank;
151 int numprocs;
152 int local_myrank;
153 int local_numprocs;
154 } phastaio_file_t;
155
156 //default: Paraview disabled
157
158 typedef struct
159 {
160 int fileID;
161 int nppf, nfields;
162 int GPid;
163 int read_field_count;
164 char * masterHeader;
165 unsigned long long **offset_table;
166 unsigned long long my_offset;
167
168 }serial_file;
169
170 serial_file *SerialFile;
171 phastaio_file_t *PhastaIOActiveFiles[MAX_PHASTA_FILES];
172 int PhastaIONextActiveIndex = 0; /* indicates next index to allocate */
173
174 //CHANGE END////////////////////////////////////////////////////////////
175
176 #define swap_char(A,B) { ucTmp = A; A = B ; B = ucTmp; }
177
178 vtkstd::map< int , char* > LastHeaderKey;
179 vtkstd::vector< FILE* > fileArray;
180 vtkstd::vector< int > byte_order;
181 vtkstd::vector< int > header_type;
182 int DataSize=0;
183 int LastHeaderNotFound = 0;
184 int Wrong_Endian = 0 ;
185 int Strict_Error = 0 ;
186 int binary_format = 0;
187
188 // the caller has the responsibility to delete the returned string
StringStripper(const char istring[])189 char* vtkPhastaReader::StringStripper( const char istring[] )
190 {
191 int length = strlen( istring );
192 char* dest = new char [ length + 1 ];
193 strcpy( dest, istring );
194 dest[ length ] = '\0';
195
196 if ( char* p = strpbrk( dest, " ") )
197 {
198 *p = '\0';
199 }
200
201 return dest;
202 }
203
cscompare(const char teststring[],const char targetstring[])204 int vtkPhastaReader::cscompare( const char teststring[],
205 const char targetstring[] )
206 {
207
208 char* s1 = const_cast<char*>(teststring);
209 char* s2 = const_cast<char*>(targetstring);
210
211 while( *s1 == ' ') { s1++; }
212 while( *s2 == ' ') { s2++; }
213 while( ( *s1 )
214 && ( *s2 )
215 && ( *s2 != '?')
216 && ( tolower( *s1 )==tolower( *s2 ) ) )
217 {
218 s1++;
219 s2++;
220 while( *s1 == ' ') { s1++; }
221 while( *s2 == ' ') { s2++; }
222 }
223 if ( !( *s1 ) || ( *s1 == '?') )
224 {
225 return 1;
226 }
227 else
228 {
229 return 0;
230 }
231 }
232
isBinary(const char iotype[])233 void vtkPhastaReader::isBinary( const char iotype[] )
234 {
235 char* fname = StringStripper( iotype );
236 if ( cscompare( fname, "binary" ) )
237 {
238 binary_format = 1;
239 }
240 else
241 {
242 binary_format = 0;
243 }
244 delete [] fname;
245
246 }
247
typeSize(const char typestring[])248 size_t vtkPhastaReader::typeSize( const char typestring[] )
249 {
250 char* ts1 = StringStripper( typestring );
251
252 if ( cscompare( "integer", ts1 ) )
253 {
254 delete [] ts1;
255 return sizeof(int);
256 }
257 else if ( cscompare( "double", ts1 ) )
258 {
259 delete [] ts1;
260 return sizeof( double );
261 }
262 else if ( cscompare( "float", ts1 ) )
263 {
264 delete [] ts1;
265 return sizeof( float );
266 }
267 else
268 {
269 delete [] ts1;
270 fprintf(stderr,"unknown type : %s\n",ts1);
271 return 0;
272 }
273 }
274
readHeader(FILE * fileObject,const char phrase[],int * params,int expect)275 int vtkPhastaReader::readHeader( FILE* fileObject,
276 const char phrase[],
277 int* params,
278 int expect )
279 {
280 char* text_header;
281 char* token;
282 char Line[1024];
283 char junk;
284 int FOUND = 0 ;
285 int real_length;
286 int skip_size, integer_value;
287 int rewind_count=0;
288
289 if( !fgets( Line, 1024, fileObject ) && feof( fileObject ) )
290 {
291 rewind( fileObject );
292 clearerr( fileObject );
293 rewind_count++;
294 fgets( Line, 1024, fileObject );
295 }
296
297 while( !FOUND && ( rewind_count < 2 ) )
298 {
299 if ( ( Line[0] != '\n' ) && ( real_length = strcspn( Line, "#" )) )
300 {
301 text_header = new char [ real_length + 1 ];
302 strncpy( text_header, Line, real_length );
303 text_header[ real_length ] =static_cast<char>(NULL);
304 token = strtok ( text_header, ":" );
305 if( cscompare( phrase , token ) )
306 {
307 FOUND = 1 ;
308 token = strtok( NULL, " ,;<>" );
309 skip_size = atoi( token );
310 int i;
311 for( i=0; i < expect && ( token = strtok( NULL," ,;<>") ); i++)
312 {
313 params[i] = atoi( token );
314 }
315 if ( i < expect )
316 {
317 fprintf(stderr,"Expected # of ints not found for: %s\n",phrase );
318 }
319 }
320 else if ( cscompare(token,"byteorder magic number") )
321 {
322 if ( binary_format )
323 {
324 fread((void*)&integer_value,sizeof(int),1,fileObject);
325 fread( &junk, sizeof(char), 1 , fileObject );
326 if ( 362436 != integer_value )
327 {
328 Wrong_Endian = 1;
329 }
330 }
331 else
332 {
333 fscanf(fileObject, "%d\n", &integer_value );
334 }
335 }
336 else
337 {
338 /* some other header, so just skip over */
339 token = strtok( NULL, " ,;<>" );
340 skip_size = atoi( token );
341 if ( binary_format)
342 {
343 fseek( fileObject, skip_size, SEEK_CUR );
344 }
345 else
346 {
347 for( int gama=0; gama < skip_size; gama++ )
348 {
349 fgets( Line, 1024, fileObject );
350 }
351 }
352 }
353 delete [] text_header;
354 }
355
356 if ( !FOUND )
357 {
358 if( !fgets( Line, 1024, fileObject ) && feof( fileObject ) )
359 {
360 rewind( fileObject );
361 clearerr( fileObject );
362 rewind_count++;
363 fgets( Line, 1024, fileObject );
364 }
365 }
366 }
367
368 if ( !FOUND )
369 {
370 fprintf(stderr, "Error: Cound not find: %s\n", phrase);
371 return 1;
372 }
373 return 0;
374 }
375
SwapArrayByteOrder_(void * array,int nbytes,int nItems)376 void vtkPhastaReader::SwapArrayByteOrder_( void* array,
377 int nbytes,
378 int nItems )
379 {
380 /* This swaps the byte order for the array of nItems each
381 of size nbytes , This will be called only locally */
382 int i,j;
383 unsigned char ucTmp;
384 unsigned char* ucDst = (unsigned char*)array;
385
386 for(i=0; i < nItems; i++)
387 {
388 for(j=0; j < (nbytes/2); j++)
389 {
390 swap_char( ucDst[j] , ucDst[(nbytes - 1) - j] );
391 }
392 ucDst += nbytes;
393 }
394 }
395
396 //CHANGE///////////////////////////////////////////////////////
397
queryphmpiio_(const char filename[],int * nfields,int * nppf)398 void vtkPhastaReader::queryphmpiio_(const char filename[],int *nfields, int *nppf)
399 {
400
401 FILE * fileHandle;
402 char* fname = StringStripper( filename );
403
404 fileHandle = fopen (fname,"rb");
405 if (fileHandle == NULL ) {
406 printf("\n File %s doesn't exist! Please check!\n",fname);
407 exit(1);
408 }
409 else
410 {
411 SerialFile =(serial_file *)calloc( 1, sizeof( serial_file) );
412
413 SerialFile->masterHeader = (char *)malloc(MasterHeaderSize);
414 fread(SerialFile->masterHeader,1,MasterHeaderSize,fileHandle);
415
416 char read_out_tag[MAX_FIELDS_NAME_LENGTH];
417 char * token;
418 int magic_number;
419 memcpy( read_out_tag,
420 SerialFile->masterHeader,
421 MAX_FIELDS_NAME_LENGTH-1 );
422
423 if ( cscompare ("MPI_IO_Tag",read_out_tag) )
424 {
425 // Test endianess ...
426 memcpy ( &magic_number,
427 SerialFile->masterHeader+sizeof("MPI_IO_Tag :"),
428 sizeof(int) );
429
430 if ( magic_number != ENDIAN_TEST_NUMBER )
431 {
432 diff_endian = 1;
433 }
434
435 char version[MAX_FIELDS_NAME_LENGTH/4];
436 int mhsize;
437
438 memcpy(version,
439 SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH/2,
440 MAX_FIELDS_NAME_LENGTH/4 - 1); //TODO: why -1?
441
442 if( cscompare ("version",version) )
443 {
444 // if there is "version" tag in the file, then it is newer format
445 // read master header size from here, otherwise use default
446 // TODO: if version is "1", we know mhsize is at 3/4 place...
447
448 token = strtok(version, ":");
449 token = strtok(NULL, " ,;<>" );
450 int iversion = atoi(token);
451
452 if( iversion == 1) {
453 memcpy( &mhsize,
454 SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH/4*3 + sizeof("mhsize : ")-1,
455 sizeof(int));
456 if ( diff_endian)
457 SwapArrayByteOrder_(&mhsize, sizeof(int), 1);
458 free(SerialFile->masterHeader);
459 SerialFile->masterHeader = (char *)malloc(mhsize);
460 fseek(fileHandle, 0, SEEK_SET);
461 fread(SerialFile->masterHeader,1,mhsize,fileHandle);
462 }
463 //TODO: check if this is a valid int??
464 MasterHeaderSize = mhsize;
465 }
466 else { // else it's version 0's format w/o version tag, implicating MHSize=4M
467 MasterHeaderSize = DefaultMHSize;
468 //printf("-----> version = 0; mhsize = %d\n", MasterHeaderSize);
469 }
470
471 // END OF CHANGE FOR VERSION
472 //
473 memcpy( read_out_tag,
474 SerialFile->masterHeader+MAX_FIELDS_NAME_LENGTH+1,
475 MAX_FIELDS_NAME_LENGTH );
476
477 // Read in # fields ...
478 token = strtok ( read_out_tag, ":" );
479 token = strtok( NULL," ,;<>" );
480 *nfields = atoi( token );
481 SerialFile->nfields=*nfields;
482
483 memcpy( read_out_tag,
484 SerialFile->masterHeader+
485 *nfields * MAX_FIELDS_NAME_LENGTH +
486 MAX_FIELDS_NAME_LENGTH * 2,
487 MAX_FIELDS_NAME_LENGTH);
488
489 token = strtok ( read_out_tag, ":" );
490 token = strtok( NULL," ,;<>" );
491 *nppf = atoi( token );
492 SerialFile->nppf=*nppf;
493 }
494 else
495 {
496 printf("The file you opened is not new format, please check!\n");
497 }
498
499 fclose(fileHandle);
500 }
501 delete [] fname;
502
503 }
504
finalizephmpiio_(int * fileDescriptor)505 void vtkPhastaReader::finalizephmpiio_( int *fileDescriptor )
506 {
507 //printf("total open time is %lf\n", opentime_total);
508 // free master header, offset table [][], and serial file struc
509 free( SerialFile->masterHeader);
510 int j;
511 for ( j = 0; j < SerialFile->nfields; j++ )
512 {
513 free( SerialFile->offset_table[j] );
514 }
515 free( SerialFile->offset_table);
516 free( SerialFile );
517 }
518
StrReverse(char * str)519 char* StrReverse(char* str)
520 {
521 char *temp, *ptr;
522 int len, i;
523
524 temp=str;
525 for(len=0; *temp !='\0';temp++, len++);
526
527 ptr=(char*)malloc(sizeof(char)*(len+1));
528
529 for(i=len-1; i>=0; i--)
530 ptr[len-i-1]=str[i];
531
532 ptr[len]='\0';
533 return ptr;
534 }
535
536
537 //CHANGE END//////////////////////////////////////////////////
538
openfile(const char filename[],const char mode[],int * fileDescriptor)539 void vtkPhastaReader::openfile( const char filename[],
540 const char mode[],
541 int* fileDescriptor )
542 //CHANGE////////////////////////////////////////////////////
543 {
544 //printf("in open(): counter = %ld\n", counter++);
545
546 FILE* file=NULL ;
547 *fileDescriptor = 0;
548 char* fname = StringStripper( filename );
549 char* imode = StringStripper( mode );
550
551 int string_length = strlen( fname );
552 char* buffer = (char*) malloc ( string_length+1 );
553 strcpy ( buffer, fname );
554 buffer[ string_length ] = '\0';
555
556 char* tempbuf = StrReverse(buffer);
557 free(buffer);
558 buffer = tempbuf;
559
560 //printf("buffer is %s\n",buffer);
561
562 char* st2 = strtok ( buffer, "." );
563 //st2 = strtok (NULL, ".");
564
565 //printf("st2 is %s\n",st2);
566
567 string_length = strlen(st2);
568 char* buffer2 = (char*)malloc(string_length+1);
569 strcpy(buffer2,st2);
570 buffer2[string_length]='\0';
571
572 char* tempbuf2 = StrReverse(buffer2);
573 free(buffer2);
574 buffer2 = tempbuf2;
575 //printf("buffer2 is %s\n",buffer2);
576
577 SerialFile->fileID = atoi(buffer2);
578 if ( char* p = strpbrk(buffer, "@") )
579 *p = '\0';
580
581 startTimer(&start);
582 if ( cscompare( "read", imode ) ) file = fopen(fname, "rb" );
583 else if( cscompare( "write", imode ) ) file = fopen(fname, "wb" );
584 else if( cscompare( "append", imode ) ) file = fopen(fname, "ab" );
585 endTimer(&end);
586 computeTime(&start, &end);
587
588
589 if ( !file ){
590 fprintf(stderr,"unable to open file : %s\n",fname ) ;
591 } else {
592 fileArray.push_back( file );
593 byte_order.push_back( false );
594 header_type.push_back( sizeof(int) );
595 *fileDescriptor = fileArray.size();
596 }
597
598 ////////////////////////////////////////////////
599 //unsigned long long **header_table;
600 SerialFile->offset_table = ( unsigned long long ** )calloc(SerialFile->nfields,
601 sizeof(unsigned long long *));
602
603 int j;
604 for ( j = 0; j < SerialFile->nfields; j++ )
605 {
606 SerialFile->offset_table[j]=( unsigned long long * ) calloc( SerialFile->nppf ,
607 sizeof( unsigned long long));
608 }
609
610 // Read in the offset table ...
611 for ( j = 0; j < SerialFile->nfields; j++ )
612 {
613
614 memcpy( SerialFile->offset_table[j],
615 SerialFile->masterHeader +
616 VERSION_INFO_HEADER_SIZE +
617 j * SerialFile->nppf * sizeof(unsigned long long),
618 SerialFile->nppf * sizeof(unsigned long long) );
619
620 if(diff_endian) {
621 SwapArrayByteOrder_( SerialFile->offset_table[j],
622 sizeof(unsigned long long int),
623 SerialFile->nppf);
624 }
625 // Swap byte order if endianess is different ...
626 /*if ( PhastaIOActiveFiles[i]->Wrong_Endian )
627 {
628 SwapArrayByteOrder_( PhastaIOActiveFiles[i]->my_read_table[j],
629 sizeof(long long int),
630 PhastaIOActiveFiles[i]->nppp );
631 }
632 */
633 }
634
635 ////////////////////////////////////////////////
636 delete [] fname;
637 delete [] imode;
638 //free(fname);
639 //free(imode);
640 free(buffer);
641 free(buffer2);
642 }
643
644 //CHANGE END////////////////////////////////////////////////
closefile(int * fileDescriptor,const char mode[])645 void vtkPhastaReader::closefile( int* fileDescriptor,
646 const char mode[] )
647 //CHANGE///////////////////////////////////////////////
648 {
649 char* imode = StringStripper( mode );
650
651 if( cscompare( "write", imode )
652 || cscompare( "append", imode ) ) {
653 fflush( fileArray[ *fileDescriptor - 1 ] );
654 }
655
656 fclose( fileArray[ *fileDescriptor - 1 ] );
657 delete [] imode;
658 }
659
660
661 //CHANGE END///////////////////////////////////////////
662
663
readheader(int * fileDescriptor,const char keyphrase[],void * valueArray,int * nItems,const char datatype[],const char iotype[])664 void vtkPhastaReader::readheader( int* fileDescriptor,
665 const char keyphrase[],
666 void* valueArray,
667 int* nItems,
668 const char datatype[],
669 const char iotype[] )
670 //CHANGE////////////////////////////////////////////////////
671 {
672 int filePtr = *fileDescriptor - 1;
673 FILE* fileObject;
674 int* valueListInt;
675
676 if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) {
677 fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor);
678 fprintf(stderr,"openfile_ function has to be called before \n") ;
679 fprintf(stderr,"acessing the file\n ") ;
680 fprintf(stderr,"fatal error: cannot continue, returning out of call\n");
681 return;
682 }
683
684 LastHeaderKey[ filePtr ] = const_cast< char* >( keyphrase );
685 LastHeaderNotFound = false;
686
687 fileObject = fileArray[ filePtr ] ;
688 Wrong_Endian = byte_order[ filePtr ];
689
690 isBinary( iotype );
691 typeSize( datatype ); //redundant call, just avoid a compiler warning.
692
693 // right now we are making the assumption that we will only write integers
694 // on the header line.
695
696 valueListInt = static_cast< int* >( valueArray );
697
698 /////////////////////////////////////////////////////////
699 int j;
700 bool FOUND = false ;
701 unsigned int skip_size;
702 char * token;
703 char readouttag[MAX_FIELDS_NUMBER][MAX_FIELDS_NAME_LENGTH];
704
705 int string_length = strlen( keyphrase );
706 char* buffer = (char*) malloc ( string_length+1 );
707 strcpy ( buffer, keyphrase );
708 buffer[ string_length ] = '\0';
709
710 char* st2 = strtok ( buffer, "@" );
711 st2 = strtok (NULL, "@");
712 SerialFile->GPid = atoi(st2);
713 if ( char* p = strpbrk(buffer, "@") )
714 *p = '\0';
715
716 //printf("field is %s and nfields is %d\n",keyphrase,SerialFile->nfields);
717
718 for ( j = 0; j<SerialFile->nfields; j++ )
719 {
720 memcpy( readouttag[j],
721 SerialFile->masterHeader + j*MAX_FIELDS_NAME_LENGTH+MAX_FIELDS_NAME_LENGTH*2+1,
722 MAX_FIELDS_NAME_LENGTH-1 );
723 }
724
725 for ( j = 0; j<SerialFile->nfields; j++ )
726 {
727 token = strtok ( readouttag[j], ":" );
728
729 if ( cscompare( buffer, token ) )
730 {
731 SerialFile->read_field_count = j;
732 FOUND = true;
733 break;
734 }
735 }
736 if (!FOUND)
737 {
738 printf("Not found %s \n",keyphrase);
739 return;
740 }
741
742 int read_part_count = SerialFile->GPid - ( SerialFile->fileID - 1 ) * SerialFile->nppf - 1;
743 SerialFile->my_offset = SerialFile->offset_table[SerialFile->read_field_count][read_part_count];
744
745
746 //printf("GP id is %d and fileID is %d and nppf is %d; ",SerialFile->GPid,SerialFile->fileID,SerialFile->nppf);
747 //printf("read field count is %d and read part count is %d; ",SerialFile->read_field_count,read_part_count);
748
749 char read_out_header[MAX_FIELDS_NAME_LENGTH];
750 fseek(fileObject, SerialFile->my_offset+1, SEEK_SET);
751 fread( read_out_header, 1, MAX_FIELDS_NAME_LENGTH-1, fileObject );
752
753
754 token = strtok ( read_out_header, ":" );
755
756 if( cscompare( keyphrase , token ) )
757 {
758 FOUND = true ;
759 token = strtok( NULL, " ,;<>" );
760 skip_size = atoi( token );
761 for( j=0; j < *nItems && ( token = strtok( NULL," ,;<>") ); j++ )
762 valueListInt[j] = atoi( token );
763 //printf("$$Keyphrase is %s Value list [0] is %d \n",keyphrase,valueListInt[0] );
764 if ( j < *nItems )
765 {
766 fprintf( stderr, "Expected # of ints not found for: %s\n", keyphrase );
767 }
768 }
769
770 /////////////////////////////////////////////////////////
771
772 byte_order[ filePtr ] = Wrong_Endian ;
773
774 //if ( ierr ) LastHeaderNotFound = true;
775
776 free(buffer);
777
778 return;
779 }
780
781 //CHANGE END////////////////////////////////////////////////
782
readdatablock(int * fileDescriptor,const char keyphrase[],void * valueArray,int * nItems,const char datatype[],const char iotype[])783 void vtkPhastaReader::readdatablock( int* fileDescriptor,
784 const char keyphrase[],
785 void* valueArray,
786 int* nItems,
787 const char datatype[],
788 const char iotype[] )
789 //CHANGE//////////////////////////////////////////////////////
790 {
791
792 int filePtr = *fileDescriptor - 1;
793 FILE* fileObject;
794 char junk;
795
796 if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) {
797 fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor);
798 fprintf(stderr,"openfile_ function has to be called before \n") ;
799 fprintf(stderr,"acessing the file\n ") ;
800 fprintf(stderr,"fatal error: cannot continue, returning out of call\n");
801 return;
802 }
803
804 // error check..
805 // since we require that a consistant header always preceed the data block
806 // let us check to see that it is actually the case.
807
808 if ( ! cscompare( LastHeaderKey[ filePtr ], keyphrase ) ) {
809 fprintf(stderr, "Header not consistant with data block\n");
810 fprintf(stderr, "Header: %s\n", LastHeaderKey[ filePtr ] );
811 fprintf(stderr, "DataBlock: %s\n ", keyphrase );
812 fprintf(stderr, "Please recheck read sequence \n");
813 if( Strict_Error ) {
814 fprintf(stderr, "fatal error: cannot continue, returning out of call\n");
815 return;
816 }
817 }
818
819 if ( LastHeaderNotFound ) return;
820
821 fileObject = fileArray[ filePtr ];
822 Wrong_Endian = byte_order[ filePtr ];
823 //printf("in readdatablock(): wrong_endian = %d\n", Wrong_Endian);
824
825 size_t type_size = typeSize( datatype );
826 int nUnits = *nItems;
827 isBinary( iotype );
828
829 if ( binary_format ) {
830 fseek(fileObject, SerialFile->my_offset+DB_HEADER_SIZE, SEEK_SET);
831
832 fread( valueArray, type_size, nUnits, fileObject );
833 //fread( &junk, sizeof(char), 1 , fileObject );
834 //if ( Wrong_Endian ) SwapArrayByteOrder_( valueArray, type_size, nUnits );
835 if ( diff_endian )
836 SwapArrayByteOrder_( valueArray, type_size, nUnits ); // fj
837 } else {
838
839 char* ts1 = StringStripper( datatype );
840 if ( cscompare( "integer", ts1 ) ) {
841 for( int n=0; n < nUnits ; n++ )
842 fscanf(fileObject, "%d\n",(int*)((int*)valueArray+n) );
843 } else if ( cscompare( "double", ts1 ) ) {
844 for( int n=0; n < nUnits ; n++ )
845 fscanf(fileObject, "%lf\n",(double*)((double*)valueArray+n) );
846 }
847 delete [] ts1;
848 }
849 return;
850 }
851
852 // End of copy from phastaIO
853
854
vtkPhastaReader()855 vtkPhastaReader::vtkPhastaReader()
856 {
857 //this->DebugOn(); // TODO: comment out this line to turn off debug
858 this->GeometryFileName = NULL;
859 this->FieldFileName = NULL;
860 this->SetNumberOfInputPorts(0);
861 this->Internal = new vtkPhastaReaderInternal;
862
863 //////////
864 this->Parser = 0;
865 this->FileName = 0;
866 //////////
867
868 }
869
~vtkPhastaReader()870 vtkPhastaReader::~vtkPhastaReader()
871 {
872 if (this->GeometryFileName)
873 {
874 delete [] this->GeometryFileName;
875 }
876 if (this->FieldFileName)
877 {
878 delete [] this->FieldFileName;
879 }
880 delete this->Internal;
881
882 ////////////////////
883 if (this->Parser)
884 this->Parser->Delete();
885 ////////////////////
886
887 }
888
ClearFieldInfo()889 void vtkPhastaReader::ClearFieldInfo()
890 {
891 this->Internal->FieldInfoMap.clear();
892 }
893
SetFieldInfo(const char * paraviewFieldTag,const char * phastaFieldTag,int index,int numOfComps,int dataDependency,const char * dataType)894 void vtkPhastaReader::SetFieldInfo(const char* paraviewFieldTag,
895 const char* phastaFieldTag,
896 int index,
897 int numOfComps,
898 int dataDependency,
899 const char* dataType)
900 {
901 //printf("In P setfino\n");
902
903 //CHANGE/////////////////////////
904 partID_counter=0;
905
906 //CHANGE END/////////////////////
907
908 vtkPhastaReaderInternal::FieldInfo &info =
909 this->Internal->FieldInfoMap[paraviewFieldTag];
910
911 info.PhastaFieldTag = phastaFieldTag;
912 info.StartIndexInPhastaArray = index;
913 info.NumberOfComponents = numOfComps;
914 info.DataDependency = dataDependency;
915 info.DataType = dataType;
916 }
917
RequestData(vtkInformation *,vtkInformationVector **,vtkInformationVector * outputVector)918 int vtkPhastaReader::RequestData(vtkInformation*,
919 vtkInformationVector**,
920 vtkInformationVector* outputVector)
921 {
922 vtkDebugMacro("In P RequestData");
923
924 int firstVertexNo = 0;
925 int fvn = 0;
926 int noOfNodes, noOfCells, noOfDatas;
927
928
929 // get the data object
930 //TODO Just Testing
931 vtkSmartPointer<vtkInformation> outInfo =
932 outputVector->GetInformationObject(0);
933
934 //change/////// This part not working
935 // get the current piece being requested
936 int piece =
937 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
938
939 //printf("piece is %d\n",piece);
940
941 //partID_counter++;
942 partID_counter=PART_ID;
943 //change end////////////////
944
945 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
946 outInfo->Get(vtkDataObject::DATA_OBJECT()));
947
948 //printf("This geom file is %s\n",this->GeometryFileName);
949
950
951 int numPieces=NUM_PIECES, numFiles=NUM_FILES, timeStep=TIME_STEP;
952 //int numPieces=8, numFiles=4, timeStep=50400;
953
954 int numPiecesPerFile = numPieces/numFiles;
955 int fileID;
956 fileID = int((partID_counter-1)/numPiecesPerFile)+1;
957
958 vtkDebugMacro(<< "FILE_PATH: " << FILE_PATH);
959 // FILE_PATH is set to be the path of .pht file ?
960 //sprintf(this->FieldFileName,"%s%s.%d.%d",FILE_PATH,"/restart-dat",timeStep,fileID); // this is Ning's version
961
962 // the file id of fieldfilename need to be changed to file id now
963 char* str = this->FieldFileName;
964 for (int i = strlen(this->FieldFileName); i >= 0; i--) {
965 if(str[i] != '.')
966 str[i] = 0;
967 else
968 break;
969 }
970 sprintf(str, "%s%d", str, FILE_ID);
971 vtkDebugMacro(<<"tweaked FieldFileName="<<this->FieldFileName);
972 ///////////////////////////////////////////////////////////
973
974 vtkPoints *points;
975
976 output->Allocate(10000, 2100);
977
978 points = vtkPoints::New();
979
980 vtkDebugMacro(<<"Reading Phasta file...");
981
982 if(!this->GeometryFileName || !this->FieldFileName )
983 {
984 vtkErrorMacro(<<"All input parameters not set.");
985 return 0;
986 }
987 vtkDebugMacro(<< "Updating ensa with ....");
988 vtkDebugMacro(<< "Geom File : " << this->GeometryFileName);
989 vtkDebugMacro(<< "Field File : " << this->FieldFileName);
990
991 fvn = firstVertexNo;
992 this->ReadGeomFile(this->GeometryFileName, firstVertexNo, points, noOfNodes, noOfCells);
993 /* set the points over here, this is because vtkUnStructuredGrid
994 only insert points once, next insertion overwrites the previous one */
995 // acbauer is not sure why the above comment is about...
996 output->SetPoints(points);
997 points->Delete();
998
999 if (!this->Internal->FieldInfoMap.size())
1000 {
1001 vtkDataSetAttributes* field = output->GetPointData();
1002 this->ReadFieldFile(this->FieldFileName, fvn, field, noOfNodes);
1003 }
1004 else
1005 {
1006 this->ReadFieldFile(this->FieldFileName, fvn, output, noOfDatas);
1007 }
1008
1009 // if there exists point arrays called coordsX, coordsY and coordsZ,
1010 // create another array of point data and set the output to use this
1011 vtkPointData* pointData = output->GetPointData();
1012 vtkDoubleArray* coordsX = vtkDoubleArray::SafeDownCast(
1013 pointData->GetArray("coordsX"));
1014 vtkDoubleArray* coordsY = vtkDoubleArray::SafeDownCast(
1015 pointData->GetArray("coordsY"));
1016 vtkDoubleArray* coordsZ = vtkDoubleArray::SafeDownCast(
1017 pointData->GetArray("coordsZ"));
1018 if(coordsX && coordsY && coordsZ)
1019 {
1020 vtkIdType numPoints = output->GetPoints()->GetNumberOfPoints();
1021 if(numPoints != coordsX->GetNumberOfTuples() ||
1022 numPoints != coordsY->GetNumberOfTuples() ||
1023 numPoints != coordsZ->GetNumberOfTuples() )
1024 {
1025 vtkWarningMacro("Wrong number of points for moving mesh. Using original points.");
1026 return 0;
1027 }
1028 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
1029 points->DeepCopy(output->GetPoints());
1030 for(vtkIdType i=0;i<numPoints;i++)
1031 {
1032 points->SetPoint(i, coordsX->GetValue(i), coordsY->GetValue(i),
1033 coordsZ->GetValue(i));
1034 }
1035 output->SetPoints(points);
1036 }
1037
1038 //printf("end of RequestData():");
1039 //printf("counter = %ld\n", counter++);
1040
1041 vtkDebugMacro("end of P RequestData() # of pieces is "<<numPieces << ", partID_counter = "<< partID_counter);
1042 vtkDebugMacro("-------> total open time is "<< opentime_total);
1043
1044 return 1;
1045 }
1046
1047 /* firstVertexNo is useful when reading multiple geom files and coalescing
1048 them into one, ReadGeomfile can then be called repeatedly from Execute with
1049 firstVertexNo forming consecutive series of vertex numbers */
1050
ReadGeomFile(char * geomFileName,int & firstVertexNo,vtkPoints * points,int & num_nodes,int & num_cells)1051 void vtkPhastaReader::ReadGeomFile(char* geomFileName,
1052 int &firstVertexNo,
1053 vtkPoints *points,
1054 int &num_nodes,
1055 int &num_cells)
1056 {
1057 vtkDebugMacro("in P ReadGeomFile(): partID="<<partID_counter);
1058
1059 /* variables for vtk */
1060 vtkUnstructuredGrid *output = this->GetOutput();
1061 double *coordinates;
1062 vtkIdType *nodes;
1063 int cell_type;
1064
1065 // int num_tpblocks;
1066 /* variables for the geom data file */
1067 /* nodal information */
1068 // int byte_order;
1069 //int data[11], data1[7];
1070 int dim;
1071 int num_int_blocks;
1072 double *pos;
1073 //int *nlworkdata;
1074 /* element information */
1075 int num_elems,num_vertices,num_per_line;
1076 int *connectivity = NULL;
1077
1078
1079 /* misc variables*/
1080 int i, j,k,item;
1081 int geomfile;
1082
1083 //CHANGE///////////////////////////////////////
1084 int getNfiles, getNPPF;
1085 queryphmpiio_(geomFileName,&getNfiles,&getNPPF);
1086 char fieldName[255];
1087
1088
1089 //CHANGE END///////////////////////////////////
1090
1091 openfile(geomFileName,"read",&geomfile);
1092 //geomfile = fopen(GeometryFileName,"rb");
1093
1094 if(!geomfile)
1095 {
1096 vtkErrorMacro(<<"Cannot open file " << geomFileName);
1097 //return;
1098 }
1099
1100 int expect;
1101 int array[10];
1102 expect = 1;
1103
1104 /* read number of nodes */
1105
1106 ///CHANGE/////////////////////////////////////////////////////
1107 bzero((void*)fieldName,255);
1108 sprintf(fieldName,"%s@%d","number of nodes",partID_counter);
1109 ///CHANGE END//////////////////////////////////////////////////
1110 readheader(&geomfile,fieldName,array,&expect,"integer","binary");
1111 //readheader(&geomfile,"number of nodes",array,&expect,"integer","binary");
1112 vtkDebugMacro("after readheader(), fieldName=" << fieldName << ", geomfile (file desc) = " << geomfile);
1113 num_nodes = array[0];
1114
1115 /* read number of elements */
1116
1117 ///CHANGE/////////////////////////////////////////////////////
1118 bzero((void*)fieldName,255);
1119 sprintf(fieldName,"%s@%d","number of interior elements",partID_counter);
1120 ///CHANGE END//////////////////////////////////////////////////
1121 readheader(&geomfile,fieldName,array,&expect,"integer","binary");
1122
1123 /*readheader(&geomfile,
1124 "number of interior elements",
1125 array,
1126 &expect,
1127 "integer",
1128 "binary");
1129 */
1130 num_elems = array[0];
1131 num_cells = array[0];
1132
1133 /* read number of interior */
1134
1135 ///CHANGE/////////////////////////////////////////////////////
1136 bzero((void*)fieldName,255);
1137 sprintf(fieldName,"%s@%d","number of interior tpblocks",partID_counter);
1138 ///CHANGE END//////////////////////////////////////////////////
1139 readheader(&geomfile,fieldName,array,&expect,"integer","binary");
1140
1141
1142 /*readheader(&geomfile,
1143 "number of interior tpblocks",
1144 array,
1145 &expect,
1146 "integer",
1147 "binary");
1148 */
1149 num_int_blocks = array[0];
1150
1151 vtkDebugMacro ( << "Nodes: " << num_nodes
1152 << "Elements: " << num_elems
1153 << "tpblocks: " << num_int_blocks );
1154
1155 /* read coordinates */
1156 expect = 2;
1157
1158 ///CHANGE/////////////////////////////////////////////////////
1159 bzero((void*)fieldName,255);
1160 sprintf(fieldName,"%s@%d","co-ordinates",partID_counter);
1161 ///CHANGE END//////////////////////////////////////////////////
1162 readheader(&geomfile,fieldName,array,&expect,"double","binary");
1163
1164
1165 //readheader(&geomfile,"co-ordinates",array,&expect,"double","binary");
1166 // TEST *******************
1167 num_nodes=array[0];
1168 // TEST *******************
1169 if(num_nodes !=array[0])
1170 {
1171 vtkErrorMacro(<<"Ambigous information in geom.data file, number of nodes does not match the co-ordinates size. Nodes: " << num_nodes << " Coordinates: " << array[0]);
1172 return;
1173 }
1174 dim = array[1];
1175
1176
1177 /* read the coordinates */
1178
1179 coordinates = new double [dim];
1180 if(coordinates == NULL)
1181 {
1182 vtkErrorMacro(<<"Unable to allocate memory for nodal info");
1183 return;
1184 }
1185
1186 pos = new double [num_nodes*dim];
1187 if(pos == NULL)
1188 {
1189 vtkErrorMacro(<<"Unable to allocate memory for nodal info");
1190 return;
1191 }
1192
1193 item = num_nodes*dim;
1194
1195 //CHANGE
1196 //readdatablock(&geomfile,"co-ordinates",pos,&item,"double","binary");
1197 //CHANGE END
1198 readdatablock(&geomfile,fieldName,pos,&item,"double","binary");
1199
1200 for(i=0;i<num_nodes;i++)
1201 {
1202 for(j=0;j<dim;j++)
1203 {
1204 coordinates[j] = pos[j*num_nodes + i];
1205 }
1206 switch(dim)
1207 {
1208 case 1:
1209 points->InsertPoint(i+firstVertexNo,coordinates[0],0,0);
1210 break;
1211 case 2:
1212 points->InsertPoint(i+firstVertexNo,coordinates[0],coordinates[1],0);
1213 break;
1214 case 3:
1215 points->InsertNextPoint(coordinates);
1216 break;
1217 default:
1218 vtkErrorMacro(<<"Unrecognized dimension in "<< geomFileName)
1219 return;
1220 }
1221 }
1222
1223 /* read the connectivity information */
1224 expect = 7;
1225
1226 for(k=0;k<num_int_blocks;k++)
1227 {
1228
1229 ///CHANGE/////////////////////////////////////////////////////
1230 bzero((void*)fieldName,255);
1231 sprintf(fieldName,"%s%d@%d","connectivity interior",k+1,partID_counter);
1232 ///CHANGE END//////////////////////////////////////////////////
1233
1234 readheader(&geomfile,
1235 fieldName,
1236 array,
1237 &expect,
1238 "integer",
1239 "binary");
1240
1241 /*readheader(&geomfile,
1242 "connectivity interior",
1243 array,
1244 &expect,
1245 "integer",
1246 "binary");
1247 */
1248
1249 /* read information about the block*/
1250 num_elems = array[0];
1251 num_vertices = array[1];
1252 num_per_line = array[3];
1253 connectivity = new int [num_elems*num_per_line];
1254
1255 if(connectivity == NULL)
1256 {
1257 vtkErrorMacro(<<"Unable to allocate memory for connectivity info");
1258 return;
1259 }
1260
1261 item = num_elems*num_per_line;
1262 /*readdatablock(&geomfile,
1263 "connectivity interior",
1264 connectivity,
1265 &item,
1266 "integer",
1267 "binary");
1268 */
1269
1270 readdatablock(&geomfile,
1271 fieldName,
1272 connectivity,
1273 &item,
1274 "integer",
1275 "binary");
1276
1277 /* insert cells */
1278 for(i=0;i<num_elems;i++)
1279 {
1280 nodes = new vtkIdType[num_vertices];
1281
1282 //connectivity starts from 1 so node[j] will never be -ve
1283 for(j=0;j<num_vertices;j++)
1284 {
1285 nodes[j] = connectivity[i+num_elems*j] + firstVertexNo - 1;
1286 }
1287
1288 /* 1 is subtracted from the connectivity info to reflect that in vtk
1289 vertex numbering start from 0 as opposed to 1 in geomfile */
1290
1291 // find out element type
1292 switch(num_vertices)
1293 {
1294 case 4:
1295 cell_type = VTK_TETRA;
1296 break;
1297 case 5:
1298 cell_type = VTK_PYRAMID;
1299 break;
1300 case 6:
1301 cell_type = VTK_WEDGE;
1302 break;
1303 case 8:
1304 cell_type = VTK_HEXAHEDRON;
1305
1306 break;
1307 default:
1308 vtkErrorMacro(<<"Unrecognized CELL_TYPE in "<< geomFileName)
1309 return;
1310 }
1311
1312 /* insert the element */
1313 output->InsertNextCell(cell_type,num_vertices,nodes);
1314 delete [] nodes;
1315 }
1316 }
1317 // update the firstVertexNo so that next slice/partition can be read
1318 firstVertexNo = firstVertexNo + num_nodes;
1319
1320 // clean up
1321 closefile(&geomfile,"read");
1322 finalizephmpiio_(&geomfile);
1323 delete [] coordinates;
1324 delete [] pos;
1325 delete [] connectivity;
1326
1327 } // end of ReadGeomFile
1328
ReadFieldFile(char * fieldFileName,int,vtkDataSetAttributes * field,int & noOfNodes)1329 void vtkPhastaReader::ReadFieldFile(char* fieldFileName,
1330 int,
1331 vtkDataSetAttributes *field,
1332 int &noOfNodes)
1333 {
1334
1335 vtkDebugMacro("In P ReadFieldFile (vtkDataSetAttr), readheader etc calls need to be updated..");
1336
1337 int i, j;
1338 int item;
1339 double *data;
1340 int fieldfile;
1341
1342 int getNfiles, getNPPF;
1343 queryphmpiio_(fieldFileName,&getNfiles,&getNPPF);
1344
1345 openfile(fieldFileName,"read",&fieldfile);
1346 //fieldfile = fopen(FieldFileName,"rb");
1347
1348 if(!fieldfile)
1349 {
1350 vtkErrorMacro(<<"Cannot open file " << FieldFileName)
1351 return;
1352 }
1353 int array[10], expect;
1354
1355 /* read the solution */
1356 vtkDoubleArray* pressure = vtkDoubleArray::New();
1357 pressure->SetName("pressure");
1358 vtkDoubleArray* velocity = vtkDoubleArray::New();
1359 velocity->SetName("velocity");
1360 velocity->SetNumberOfComponents(3);
1361 vtkDoubleArray* temperature = vtkDoubleArray::New();
1362 temperature->SetName("temperature");
1363
1364 expect = 3;
1365 readheader(&fieldfile,"solution",array,&expect,"double","binary");
1366 noOfNodes = array[0];
1367 this->NumberOfVariables = array[1];
1368
1369 vtkDoubleArray* sArrays[4];
1370 for (i=0; i<4; i++)
1371 {
1372 sArrays[i] = 0;
1373 }
1374 item = noOfNodes*this->NumberOfVariables;
1375 data = new double[item];
1376 if(data == NULL)
1377 {
1378 vtkErrorMacro(<<"Unable to allocate memory for field info");
1379 return;
1380 }
1381
1382 readdatablock(&fieldfile,"solution",data,&item,"double","binary");
1383
1384 for (i=5; i<this->NumberOfVariables; i++)
1385 {
1386 int idx=i-5;
1387 sArrays[idx] = vtkDoubleArray::New();
1388 vtksys_ios::ostringstream aName;
1389 aName << "s" << idx+1 << ends;
1390 sArrays[idx]->SetName(aName.str().c_str());
1391 sArrays[idx]->SetNumberOfTuples(noOfNodes);
1392 }
1393
1394 pressure->SetNumberOfTuples(noOfNodes);
1395 velocity->SetNumberOfTuples(noOfNodes);
1396 temperature->SetNumberOfTuples(noOfNodes);
1397 for(i=0;i<noOfNodes;i++)
1398 {
1399 pressure->SetTuple1(i, data[i]);
1400 velocity->SetTuple3(i,
1401 data[noOfNodes + i],
1402 data[2*noOfNodes + i],
1403 data[3*noOfNodes + i]);
1404 temperature->SetTuple1(i, data[4*noOfNodes + i]);
1405 for(j=5; j<this->NumberOfVariables; j++)
1406 {
1407 sArrays[j-5]->SetTuple1(i, data[j*noOfNodes + i]);
1408 }
1409 }
1410
1411 field->AddArray(pressure);
1412 field->SetActiveScalars("pressure");
1413 pressure->Delete();
1414 field->AddArray(velocity);
1415 field->SetActiveVectors("velocity");
1416 velocity->Delete();
1417 field->AddArray(temperature);
1418 temperature->Delete();
1419
1420 for (i=5; i<this->NumberOfVariables; i++)
1421 {
1422 int idx=i-5;
1423 field->AddArray(sArrays[idx]);
1424 sArrays[idx]->Delete();
1425 }
1426
1427 // clean up
1428 closefile(&fieldfile,"read");
1429 finalizephmpiio_(&fieldfile);
1430 delete [] data;
1431
1432 } //closes ReadFieldFile (vtkDataSetAttr)
1433
1434
ReadFieldFile(char * fieldFileName,int,vtkUnstructuredGrid * output,int & noOfDatas)1435 void vtkPhastaReader::ReadFieldFile(char* fieldFileName,
1436 int,
1437 vtkUnstructuredGrid *output,
1438 int &noOfDatas)
1439 {
1440
1441 vtkDebugMacro("In P ReadFieldFile (vtkUnstructuredGrid): fieldFileName="<<fieldFileName << ", partID_counter=" << partID_counter);
1442
1443 int i, j, numOfVars;
1444 int item;
1445 int fieldfile;
1446
1447 //printf("In P readfieldfile 2 and %d\n",partID_counter);
1448 //printf("read filed file name is %s\n",fieldFileName);
1449
1450 //CHANGE///////////////////////////////////////
1451 int getNfiles, getNPPF;
1452 queryphmpiio_(fieldFileName,&getNfiles,&getNPPF);
1453 char fieldName[255];
1454
1455 //CHANGE END///////////////////////////////////
1456
1457 openfile(fieldFileName,"read",&fieldfile);
1458 //fieldfile = fopen(FieldFileName,"rb");
1459
1460 //printf("open handle is %d\n",fieldfile);
1461
1462 if(!fieldfile)
1463 {
1464 vtkErrorMacro(<<"Cannot open file " << FieldFileName)
1465 //return;
1466 }
1467 int array[10], expect;
1468
1469 int activeScalars = 0, activeTensors = 0;
1470
1471 vtkPhastaReaderInternal::FieldInfoMapType::iterator it = this->Internal->FieldInfoMap.begin();
1472 vtkPhastaReaderInternal::FieldInfoMapType::iterator itend = this->Internal->FieldInfoMap.end();
1473 for(; it!=itend; it++)
1474 {
1475 const char* paraviewFieldTag = it->first.c_str();
1476 const char* phastaFieldTag = it->second.PhastaFieldTag.c_str();
1477 int index = it->second.StartIndexInPhastaArray;
1478 int numOfComps = it->second.NumberOfComponents;
1479 int dataDependency = it->second.DataDependency;
1480 const char* dataType = it->second.DataType.c_str();
1481
1482
1483 vtkDataSetAttributes* field;
1484 if(dataDependency)
1485 field = output->GetCellData();
1486 else
1487 field = output->GetPointData();
1488
1489 // void *data;
1490 int dtype; // (0=double, 1=float)
1491 vtkDataArray *dataArray;
1492 /* read the field data */
1493 if(strcmp(dataType,"double")==0)
1494 {
1495 dataArray = vtkDoubleArray::New();
1496 dtype = 0;
1497 }
1498 else if(strcmp(dataType,"float")==0)
1499 {
1500 dataArray = vtkFloatArray::New();
1501 dtype=1;
1502 }
1503 else
1504 {
1505 vtkErrorMacro("Data type [" << dataType <<"] NOT supported");
1506 continue;
1507 }
1508
1509 dataArray->SetName(paraviewFieldTag);
1510 dataArray->SetNumberOfComponents(numOfComps);
1511
1512 expect = 3;
1513
1514 //printf("version 1: hello 3 in P\n");
1515
1516 ///CHANGE/////////////////////////////////////////////////////
1517 bzero((void*)fieldName,255);
1518 sprintf(fieldName,"%s@%d",phastaFieldTag,partID_counter);
1519 ///CHANGE END//////////////////////////////////////////////////
1520 readheader(&fieldfile,fieldName,array,&expect,dataType,"binary");
1521
1522 //printf("fieldfile is %d expect is %d and datatype is %s\n",fieldfile,expect,dataType);
1523
1524 //printf("original fieldname is %s and fieldname is %s\n",phastaFieldTag,fieldName);
1525 //printf("array 0 is %d and array 1 is %d\n",array[0],array[1]);
1526
1527 //readheader(&fieldfile,phastaFieldTag,array,&expect,dataType,"binary");
1528 noOfDatas = array[0];
1529 this->NumberOfVariables = array[1];
1530 numOfVars = array[1];
1531 dataArray->SetNumberOfTuples(noOfDatas);
1532
1533 //printf("hello 4 in P\n");
1534
1535 if(index<0 || index>numOfVars-1)
1536 {
1537 vtkErrorMacro("index ["<<index<<"] is out of range [num. of vars.:"<<numOfVars<<"] for field [paraview field tag:"<<paraviewFieldTag<<", phasta field tag:"<<phastaFieldTag<<"]");
1538
1539 dataArray->Delete();
1540 continue;
1541 }
1542
1543 //printf("hello 5 in P\n");
1544
1545 if(numOfComps<0 || index+numOfComps>numOfVars)
1546 {
1547 vtkErrorMacro("index ["<<index<<"] with num. of comps. ["<<numOfComps<<"] is out of range [num. of vars.:"<<numOfVars<<"] for field [paraview field tag:"<<paraviewFieldTag<<", phasta field tag:"<<phastaFieldTag<<"]");
1548
1549 dataArray->Delete();
1550 continue;
1551 }
1552
1553 //printf("hello 6 in P\n");
1554
1555 item = numOfVars*noOfDatas;
1556
1557 //printf("item is %d and dtype is %d and numofvars is %d num of data is %d\n",item,dtype,numOfVars,noOfDatas);
1558
1559 if (dtype==0)
1560 { //data is type double
1561
1562 double *data;
1563 data = new double[item];
1564
1565 if(data == NULL)
1566 {
1567 vtkErrorMacro(<<"Unable to allocate memory for field info");
1568
1569 dataArray->Delete();
1570 continue;
1571 }
1572
1573 //printf("hello 8 in P\n");
1574
1575 //////////CHANGE/////////////////////
1576 readdatablock(&fieldfile,fieldName,data,&item,dataType,"binary");
1577 ///////////CHANGE END///////////////
1578
1579 //readdatablock(&fieldfile,phastaFieldTag,data,&item,dataType,"binary");
1580
1581 //printf("hello 9 in P\n");
1582
1583 switch(numOfComps)
1584 {
1585 case 1 :
1586 {
1587 int offset = index*noOfDatas;
1588
1589 if(!activeScalars)
1590 field->SetActiveScalars(paraviewFieldTag);
1591 else
1592 activeScalars = 1;
1593 for(i=0;i<noOfDatas;i++)
1594 {
1595 dataArray->SetTuple1(i, data[offset+i]);
1596 }
1597 }
1598 break;
1599 case 3 :
1600 {
1601 int offset[3];
1602 for(j=0;j<3;j++)
1603 offset[j] = (index+j)*noOfDatas;
1604
1605 if(!activeScalars)
1606 field->SetActiveVectors(paraviewFieldTag);
1607 else
1608 activeScalars = 1;
1609 for(i=0;i<noOfDatas;i++)
1610 {
1611 dataArray->SetTuple3(i,
1612 data[offset[0]+i],
1613 data[offset[1]+i],
1614 data[offset[2]+i]);
1615 }
1616 }
1617 break;
1618 case 9 :
1619 {
1620 int offset[9];
1621 for(j=0;j<9;j++)
1622 offset[j] = (index+j)*noOfDatas;
1623
1624 if(!activeTensors)
1625 field->SetActiveTensors(paraviewFieldTag);
1626 else
1627 activeTensors = 1;
1628 for(i=0;i<noOfDatas;i++)
1629 {
1630 dataArray->SetTuple9(i,
1631 data[offset[0]+i],
1632 data[offset[1]+i],
1633 data[offset[2]+i],
1634 data[offset[3]+i],
1635 data[offset[4]+i],
1636 data[offset[5]+i],
1637 data[offset[6]+i],
1638 data[offset[7]+i],
1639 data[offset[8]+i]);
1640 }
1641 }
1642 break;
1643 default:
1644 vtkErrorMacro("number of components [" << numOfComps <<"] NOT supported");
1645
1646 dataArray->Delete();
1647 delete [] data;
1648 continue;
1649 }
1650
1651 // clean up
1652 delete [] data;
1653
1654 }
1655 else if(dtype==1)
1656 { // data is type float
1657
1658 float *data;
1659 data = new float[item];
1660
1661 if(data == NULL)
1662 {
1663 vtkErrorMacro(<<"Unable to allocate memory for field info");
1664
1665 dataArray->Delete();
1666 continue;
1667 }
1668
1669 //////////CHANGE/////////////////////
1670 readdatablock(&fieldfile,fieldName,data,&item,dataType,"binary");
1671 //////////CHANGE END/////////////////
1672
1673 //readdatablock(&fieldfile,phastaFieldTag,data,&item,dataType,"binary");
1674
1675 switch(numOfComps)
1676 {
1677 case 1 :
1678 {
1679 int offset = index*noOfDatas;
1680
1681 if(!activeScalars)
1682 field->SetActiveScalars(paraviewFieldTag);
1683 else
1684 activeScalars = 1;
1685 for(i=0;i<noOfDatas;i++)
1686 {
1687 //double tmpval = (double) data[offset+i];
1688 //dataArray->SetTuple1(i, tmpval);
1689 dataArray->SetTuple1(i, data[offset+i]);
1690 }
1691 }
1692 break;
1693 case 3 :
1694 {
1695 int offset[3];
1696 for(j=0;j<3;j++)
1697 offset[j] = (index+j)*noOfDatas;
1698
1699 if(!activeScalars)
1700 field->SetActiveVectors(paraviewFieldTag);
1701 else
1702 activeScalars = 1;
1703 for(i=0;i<noOfDatas;i++)
1704 {
1705 //double tmpval[3];
1706 //for(j=0;j<3;j++)
1707 // tmpval[j] = (double) data[offset[j]+i];
1708 //dataArray->SetTuple3(i,
1709 // tmpval[0], tmpval[1], tmpval[3]);
1710 dataArray->SetTuple3(i,
1711 data[offset[0]+i],
1712 data[offset[1]+i],
1713 data[offset[2]+i]);
1714 }
1715 }
1716 break;
1717 case 9 :
1718 {
1719 int offset[9];
1720 for(j=0;j<9;j++)
1721 offset[j] = (index+j)*noOfDatas;
1722
1723 if(!activeTensors)
1724 field->SetActiveTensors(paraviewFieldTag);
1725 else
1726 activeTensors = 1;
1727 for(i=0;i<noOfDatas;i++)
1728 {
1729 //double tmpval[9];
1730 //for(j=0;j<9;j++)
1731 // tmpval[j] = (double) data[offset[j]+i];
1732 //dataArray->SetTuple9(i,
1733 // tmpval[0],
1734 // tmpval[1],
1735 // tmpval[2],
1736 // tmpval[3],
1737 // tmpval[4],
1738 // tmpval[5],
1739 // tmpval[6],
1740 // tmpval[7],
1741 // tmpval[8]);
1742 dataArray->SetTuple9(i,
1743 data[offset[0]+i],
1744 data[offset[1]+i],
1745 data[offset[2]+i],
1746 data[offset[3]+i],
1747 data[offset[4]+i],
1748 data[offset[5]+i],
1749 data[offset[6]+i],
1750 data[offset[7]+i],
1751 data[offset[8]+i]);
1752 }
1753 }
1754 break;
1755 default:
1756 vtkErrorMacro("number of components [" << numOfComps <<"] NOT supported");
1757
1758 dataArray->Delete();
1759 delete [] data;
1760 continue;
1761 }
1762
1763 // clean up
1764 delete [] data;
1765
1766 }
1767 else
1768 {
1769 vtkErrorMacro("Data type [" << dataType <<"] NOT supported");
1770 continue;
1771 }
1772
1773
1774 field->AddArray(dataArray);
1775
1776 // clean up
1777 dataArray->Delete();
1778 //delete [] data;
1779 }
1780
1781 // close up
1782 closefile(&fieldfile,"read");
1783 finalizephmpiio_(&fieldfile);
1784
1785 }//closes ReadFieldFile (vtkUnstructuredGrid)
1786
PrintSelf(ostream & os,vtkIndent indent)1787 void vtkPhastaReader::PrintSelf(ostream& os, vtkIndent indent)
1788 {
1789 this->Superclass::PrintSelf(os,indent);
1790
1791 os << indent << "GeometryFileName: "
1792 << (this->GeometryFileName?this->GeometryFileName:"(none)")
1793 << endl;
1794 os << indent << "FieldFileName: "
1795 << (this->FieldFileName?this->FieldFileName:"(none)")
1796 << endl;
1797 }
1798