1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkVeraOutReader.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 // .NAME vtkVeraOutReader - File reader for VERA OUT HDF5 format.
16
17 #include "vtkVeraOutReader.h"
18
19 // vtkCommonCore
20 #include "vtkDataArray.h"
21 #include "vtkDataArraySelection.h"
22 #include "vtkDoubleArray.h"
23 #include "vtkFloatArray.h"
24 #include "vtkInformation.h"
25 #include "vtkInformationDoubleVectorKey.h"
26 #include "vtkInformationVector.h"
27 #include "vtkIntArray.h"
28 #include "vtkLongArray.h"
29 #include "vtkLongLongArray.h"
30 #include "vtkObjectFactory.h"
31 #include "vtkShortArray.h"
32 #include "vtkUnsignedCharArray.h"
33 #include "vtkUnsignedIntArray.h"
34 #include "vtkUnsignedShortArray.h"
35
36 // vtkCommonExecutionModel
37 #include "vtkStreamingDemandDrivenPipeline.h"
38
39 // vtkCommonDataModel
40 #include "vtkCellData.h"
41 #include "vtkFieldData.h"
42 #include "vtkRectilinearGrid.h"
43
44 // vtkhdf5
45 #define H5_USE_16_API
46 #include <vtk_hdf5.h>
47
48 #include <numeric>
49 #include <sstream>
50 #include <vector>
51
52 using namespace std;
53
54 #define VERA_MAX_DIMENSION 6
55 #define DATASET_NAME_MAX_SIZE 1024
56
57 //*****************************************************************************
58 class vtkVeraOutReader::Internals
59 {
60 public:
Internals(vtkObject * owner)61 Internals(vtkObject* owner)
62 {
63 this->Owner = owner;
64
65 this->FileId = -1;
66
67 this->NumberOfDimensions = 0;
68
69 this->NeedCoreProcessing = true;
70
71 this->APITCH = 20; // FIXME
72 this->NASSX = 4;
73 this->NASSY = 4;
74 this->NPIN = 17;
75 this->NAX = 4;
76 this->APITCH = 20;
77 this->SYMMETRY = 0;
78 this->NUMBER_OF_STATES = 0;
79 }
80 // --------------------------------------------------------------------------
~Internals()81 virtual ~Internals() { this->CloseFile(); }
82
83 // --------------------------------------------------------------------------
84
SetFileName(const char * filename)85 void SetFileName(const char* filename)
86 {
87 std::string newFileName(filename ? filename : "");
88 ;
89 if (newFileName != this->FileName)
90 {
91 this->FileName = filename;
92 this->CloseFile();
93
94 // Reset any cache
95 this->NUMBER_OF_STATES = 0;
96 this->NeedCoreProcessing = true;
97 this->CoreCellData.clear();
98 this->CellDataArraySelection->RemoveAllArrays();
99 }
100 }
101
102 // --------------------------------------------------------------------------
103
OpenFile()104 bool OpenFile()
105 {
106 if (this->FileId > -1)
107 {
108 // Already open, skip...
109 return true;
110 }
111
112 hid_t fileAccessPropListID = H5Pcreate(H5P_FILE_ACCESS);
113 if (fileAccessPropListID < 0)
114 {
115 vtkErrorWithObjectMacro(this->Owner, "Couldn't H5Pcreate");
116 return false;
117 }
118 herr_t err = H5Pset_fclose_degree(fileAccessPropListID, H5F_CLOSE_SEMI);
119 if (err < 0)
120 {
121 vtkErrorWithObjectMacro(this->Owner, "Couldn't set file close access");
122 return false;
123 }
124 if ((this->FileId = H5Fopen(this->FileName.c_str(), H5F_ACC_RDONLY, fileAccessPropListID)) < 0)
125 {
126 vtkErrorWithObjectMacro(
127 this->Owner, "Cannot be a VERA file (" << this->FileName.c_str() << ")");
128 return false;
129 }
130 H5Pclose(fileAccessPropListID);
131 return true;
132 }
133
134 //----------------------------------------------------------------------------
135
CloseFile()136 void CloseFile()
137 {
138 if (this->FileId > -1)
139 {
140 H5Fclose(this->FileId);
141 this->FileId = -1;
142 }
143 }
144
145 //----------------------------------------------------------------------------
146
LoadMetaData()147 void LoadMetaData()
148 {
149 if (this->FileId == -1)
150 {
151 return;
152 }
153
154 this->ReadCore();
155
156 // Register state fields
157 if (this->GetNumberOfTimeSteps())
158 {
159 // Open group
160 hid_t groupId = -1;
161 if ((groupId = H5Gopen(this->FileId, "/STATE_0001")) < 0)
162 {
163 vtkErrorWithObjectMacro(this->Owner, "Can't open Group /STATE_0001");
164 return;
165 }
166
167 H5G_info_t groupInfo;
168 int status = H5Gget_info(groupId, &groupInfo);
169 if (status < 0)
170 {
171 vtkErrorWithObjectMacro(
172 this->Owner, "Can't get group info for /STATE_0001 (status: " << status << ")");
173 H5Gclose(groupId);
174 return;
175 }
176
177 // Extract dataset names
178 std::vector<std::string> datasetNames;
179 char datasetName[DATASET_NAME_MAX_SIZE];
180 for (hsize_t idx = 0; idx < groupInfo.nlinks; idx++)
181 {
182 H5Lget_name_by_idx(groupId, ".", H5_INDEX_NAME, H5_ITER_INC, idx, datasetName,
183 DATASET_NAME_MAX_SIZE, H5P_DEFAULT);
184 datasetNames.push_back(datasetName);
185 }
186 H5Gclose(groupId);
187
188 // Start processing datasets
189 for (const std::string& dsName : datasetNames)
190 {
191 this->ReadDataSetDimensions("/STATE_0001", dsName.c_str());
192 if (this->NumberOfDimensions == 4 && this->Dimensions[0] == this->NPIN &&
193 this->Dimensions[1] == this->NPIN && this->Dimensions[2] == this->NAX &&
194 this->Dimensions[3] == this->NASS)
195 {
196 this->CellDataArraySelection->AddArray(dsName.c_str());
197 }
198 else if (this->NumberOfDimensions == 1 && this->Dimensions[0] == 1)
199 {
200 this->FieldDataArraySelection->AddArray(dsName.c_str());
201 }
202 }
203 }
204 }
205
206 //----------------------------------------------------------------------------
207
GetNumberOfTimeSteps()208 int GetNumberOfTimeSteps()
209 {
210 if (this->NUMBER_OF_STATES)
211 {
212 return this->NUMBER_OF_STATES;
213 }
214
215 if (this->FileId == -1)
216 {
217 return 0;
218 }
219
220 // ----------------------------------
221 // Find out how many state we have
222 // ----------------------------------
223 int count = 0;
224 int status = 0;
225 while (status >= 0)
226 {
227 count++;
228 std::ostringstream groupName;
229 groupName << "/STATE_" << std::setw(4) << std::setfill('0') << count;
230 H5Eset_auto(NULL, NULL);
231 status = H5Gget_objinfo(this->FileId, groupName.str().c_str(), 0, NULL);
232 }
233 // H5Eset_auto(NULL, NULL);
234 this->NUMBER_OF_STATES = count ? count - 1 : 0;
235 // ----------------------------------
236
237 return this->NUMBER_OF_STATES;
238 }
239
240 //----------------------------------------------------------------------------
241
ReadDataSetDimensions(const char * groupName,const char * datasetName)242 bool ReadDataSetDimensions(const char* groupName, const char* datasetName)
243 {
244 if (this->FileId == -1)
245 {
246 return false;
247 }
248
249 // Open group
250 hid_t groupId = -1;
251 if ((groupId = H5Gopen(this->FileId, groupName)) < 0)
252 {
253 vtkErrorWithObjectMacro(this->Owner, "Can't open group " << groupName);
254 return false;
255 }
256
257 // Open dataset
258 hid_t datasetId;
259 if ((datasetId = H5Dopen(groupId, datasetName)) < 0)
260 {
261 H5Gclose(groupId);
262 vtkErrorWithObjectMacro(this->Owner, "DataSet " << datasetName << " in group " << groupName
263 << " don't want to open.");
264 return false;
265 }
266
267 // Extract dataset dimensions
268 hid_t spaceId = H5Dget_space(datasetId);
269 H5Sget_simple_extent_dims(spaceId, this->Dimensions, nullptr);
270 this->NumberOfDimensions = H5Sget_simple_extent_ndims(spaceId);
271
272 // Close resources
273 H5Sclose(spaceId);
274 H5Dclose(datasetId);
275 H5Gclose(groupId);
276
277 return true;
278 }
279
280 //----------------------------------------------------------------------------
281
ReadDataSet(const char * groupName,const char * datasetName)282 vtkDataArray* ReadDataSet(const char* groupName, const char* datasetName)
283 {
284 if (!this->ReadDataSetDimensions(groupName, datasetName))
285 {
286 return nullptr;
287 }
288
289 vtkDataArray* arrayToReturn = nullptr;
290
291 // Compute total size
292 vtkIdType nbTuples = 1;
293 for (hsize_t idx = 0; idx < this->NumberOfDimensions; idx++)
294 {
295 nbTuples *= this->Dimensions[idx];
296 }
297
298 // Open group
299 hid_t groupId = -1;
300 if ((groupId = H5Gopen(this->FileId, groupName)) < 0)
301 {
302 vtkErrorWithObjectMacro(this->Owner, "Can't open group " << groupName);
303 return nullptr;
304 }
305
306 // Open dataset
307 hid_t datasetId;
308 if ((datasetId = H5Dopen(groupId, datasetName)) < 0)
309 {
310 vtkErrorWithObjectMacro(this->Owner, "DataSet " << datasetName << " in group " << groupName
311 << " don't want to open.");
312 H5Gclose(groupId);
313 return nullptr;
314 }
315
316 // Extract data type
317 hid_t tRawType = H5Dget_type(datasetId);
318 hid_t dataType = H5Tget_native_type(tRawType, H5T_DIR_ASCEND);
319 if (H5Tequal(dataType, H5T_NATIVE_FLOAT))
320 {
321 arrayToReturn = vtkFloatArray::New();
322 arrayToReturn->SetNumberOfTuples(nbTuples);
323 float* arrayPtr =
324 static_cast<float*>(vtkArrayDownCast<vtkFloatArray>(arrayToReturn)->GetPointer(0));
325 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
326 arrayPtr = nullptr;
327 }
328 else if (H5Tequal(dataType, H5T_NATIVE_DOUBLE))
329 {
330 arrayToReturn = vtkDoubleArray::New();
331 arrayToReturn->SetNumberOfTuples(nbTuples);
332 double* arrayPtr =
333 static_cast<double*>(vtkArrayDownCast<vtkDoubleArray>(arrayToReturn)->GetPointer(0));
334 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
335 arrayPtr = nullptr;
336 }
337 else if (H5Tequal(dataType, H5T_NATIVE_INT))
338 {
339 arrayToReturn = vtkIntArray::New();
340 arrayToReturn->SetNumberOfTuples(nbTuples);
341 int* arrayPtr =
342 static_cast<int*>(vtkArrayDownCast<vtkIntArray>(arrayToReturn)->GetPointer(0));
343 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
344 arrayPtr = nullptr;
345 }
346 else if (H5Tequal(dataType, H5T_NATIVE_UINT))
347 {
348 arrayToReturn = vtkUnsignedIntArray::New();
349 arrayToReturn->SetNumberOfTuples(nbTuples);
350 unsigned int* arrayPtr = static_cast<unsigned int*>(
351 vtkArrayDownCast<vtkUnsignedIntArray>(arrayToReturn)->GetPointer(0));
352 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
353 arrayPtr = nullptr;
354 }
355 else if (H5Tequal(dataType, H5T_NATIVE_SHORT))
356 {
357 arrayToReturn = vtkShortArray::New();
358 arrayToReturn->SetNumberOfTuples(nbTuples);
359 short* arrayPtr =
360 static_cast<short*>(vtkArrayDownCast<vtkShortArray>(arrayToReturn)->GetPointer(0));
361 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
362 arrayPtr = nullptr;
363 }
364 else if (H5Tequal(dataType, H5T_NATIVE_USHORT))
365 {
366 arrayToReturn = vtkUnsignedShortArray::New();
367 arrayToReturn->SetNumberOfTuples(nbTuples);
368 unsigned short* arrayPtr = static_cast<unsigned short*>(
369 vtkArrayDownCast<vtkUnsignedShortArray>(arrayToReturn)->GetPointer(0));
370 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
371 arrayPtr = nullptr;
372 }
373 else if (H5Tequal(dataType, H5T_NATIVE_UCHAR))
374 {
375
376 arrayToReturn = vtkUnsignedCharArray::New();
377 arrayToReturn->SetNumberOfTuples(nbTuples);
378 unsigned char* arrayPtr = static_cast<unsigned char*>(
379 vtkArrayDownCast<vtkUnsignedCharArray>(arrayToReturn)->GetPointer(0));
380 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
381 arrayPtr = nullptr;
382 }
383 else if (H5Tequal(dataType, H5T_NATIVE_LONG))
384 {
385 arrayToReturn = vtkLongArray::New();
386 arrayToReturn->SetNumberOfTuples(nbTuples);
387 long* arrayPtr =
388 static_cast<long*>(vtkArrayDownCast<vtkLongArray>(arrayToReturn)->GetPointer(0));
389 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
390 arrayPtr = nullptr;
391 }
392 else if (H5Tequal(dataType, H5T_NATIVE_LLONG))
393 {
394 arrayToReturn = vtkLongLongArray::New();
395 arrayToReturn->SetNumberOfTuples(nbTuples);
396 long long* arrayPtr =
397 static_cast<long long*>(vtkArrayDownCast<vtkLongLongArray>(arrayToReturn)->GetPointer(0));
398 H5Dread(datasetId, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, arrayPtr);
399 arrayPtr = nullptr;
400 }
401 else
402 {
403 vtkErrorWithObjectMacro(this->Owner, "Unknown HDF5 data type --- it is not FLOAT, "
404 << "DOUBLE, INT, UNSIGNED INT, SHORT, UNSIGNED SHORT, "
405 << "UNSIGNED CHAR, LONG, or LONG LONG.");
406 }
407
408 arrayToReturn->SetName(datasetName);
409
410 // Close what we opened
411 H5Tclose(dataType);
412 // H5Tclose( tRawType );
413 H5Dclose(datasetId);
414 H5Gclose(groupId);
415
416 return arrayToReturn;
417 }
418
419 //----------------------------------------------------------------------------
420
CreatePinFieldArray(vtkDataArray * dataSource)421 vtkDataArray* CreatePinFieldArray(vtkDataArray* dataSource)
422 {
423 if (this->CoreMap == nullptr || dataSource == nullptr)
424 {
425 return nullptr;
426 }
427
428 vtkDataArray* outputField = dataSource->NewInstance();
429 outputField->SetNumberOfTuples(this->NASSX * this->NPIN * this->NASSY * this->NPIN * this->NAX);
430 vtkIdType srcIndex = 0;
431
432 for (hsize_t sj = 0; sj < this->NASSY; sj++)
433 {
434 for (hsize_t si = 0; si < this->NASSX; si++)
435 {
436 vtkIdType assembyId = this->CoreMap->GetTuple1(si * this->NASSX + sj) - 1;
437 vtkIdType dstOffset = si * this->NPIN + sj * this->NASSX * this->NPIN * this->NPIN;
438 for (hsize_t dk = 0; dk < this->NAX; dk++)
439 {
440 for (hsize_t dj = 0; dj < this->NPIN; dj++)
441 {
442 for (hsize_t di = 0; di < this->NPIN; di++)
443 {
444 if (assembyId < 0)
445 {
446 outputField->SetTuple1(dstOffset + di + (dj * this->NASSX * this->NPIN) +
447 (dk * this->NASSX * this->NPIN * this->NASSY * this->NPIN),
448 0);
449 }
450 else
451 {
452 // Fortran ordering
453 srcIndex = assembyId + dk * this->NASS + dj * this->NASS * this->NAX +
454 di * this->NASS * this->NAX * this->NPIN;
455 if (SYMMETRY == 4)
456 {
457 srcIndex = assembyId + dk * this->NASS;
458 if (2 * si > this->NASSX)
459 {
460 srcIndex += di * this->NASS * this->NAX * this->NPIN;
461 }
462 else
463 {
464 srcIndex += (this->NPIN - di - 1) * this->NASS * this->NAX * this->NPIN;
465 }
466 if (2 * sj > this->NASSY)
467 {
468 srcIndex += dj * this->NASS * this->NAX;
469 }
470 else
471 {
472 srcIndex += (this->NPIN - dj - 1) * this->NASS * this->NAX;
473 }
474 }
475 outputField->SetTuple1(dstOffset + di + (dj * this->NASSX * this->NPIN) +
476 (dk * this->NASSX * this->NPIN * this->NASSY * this->NPIN),
477 dataSource->GetTuple1(srcIndex));
478 }
479 }
480 }
481 }
482 }
483 }
484 return outputField;
485 }
486
487 //----------------------------------------------------------------------------
488
AddDataSetNamesWithDimension(const char * groupName,hsize_t dimension,std::vector<std::string> & names)489 void AddDataSetNamesWithDimension(
490 const char* groupName, hsize_t dimension, std::vector<std::string>& names)
491 {
492 // Open group
493 hid_t groupId = -1;
494 if ((groupId = H5Gopen(this->FileId, groupName)) < 0)
495 {
496 vtkErrorWithObjectMacro(this->Owner, "Can't open group " << groupName);
497 return;
498 }
499
500 H5G_info_t groupInfo;
501 int status = H5Gget_info(groupId, &groupInfo);
502 if (status < 0)
503 {
504 vtkErrorWithObjectMacro(this->Owner, "Can't get group info for " << groupName);
505 H5Gclose(groupId);
506 return;
507 }
508
509 // Extract dataset names
510 std::vector<std::string> datasetNames;
511 char datasetName[DATASET_NAME_MAX_SIZE];
512 for (hsize_t idx = 0; idx < groupInfo.nlinks; idx++)
513 {
514 H5Lget_name_by_idx(groupId, ".", H5_INDEX_NAME, H5_ITER_INC, idx, datasetName,
515 DATASET_NAME_MAX_SIZE, H5P_DEFAULT);
516 datasetNames.push_back(datasetName);
517 }
518
519 // Start processing datasets
520 for (const std::string& dsName : datasetNames)
521 {
522 // Open dataset
523 hid_t datasetId;
524 if ((datasetId = H5Dopen(groupId, dsName.c_str())) < 0)
525 {
526 vtkErrorWithObjectMacro(this->Owner, "DataSet " << dsName.c_str() << " in group "
527 << groupName << " don't want to open.");
528 continue;
529 }
530
531 // Extract dataset dimensions
532 hid_t spaceId = H5Dget_space(datasetId);
533 H5Sget_simple_extent_dims(spaceId, this->Dimensions, nullptr);
534 this->NumberOfDimensions = H5Sget_simple_extent_ndims(spaceId);
535
536 // Add name if condition is valid
537 if (this->NumberOfDimensions == dimension)
538 {
539 names.push_back(dsName);
540 }
541 H5Sclose(spaceId);
542 H5Dclose(datasetId);
543 }
544 H5Gclose(groupId);
545 }
546
547 //----------------------------------------------------------------------------
548
ReadCore()549 void ReadCore()
550 {
551 if (!this->NeedCoreProcessing)
552 {
553 return;
554 }
555
556 // Guard further reading if file name does not change..
557 this->NeedCoreProcessing = false;
558 this->CoreCellData.clear();
559
560 // Local vars
561 vtkDataArray* dataSource;
562 vtkDataArray* outputCellArray;
563
564 // --------------------------------------------------------------------------
565 // Global variables
566 // --------------------------------------------------------------------------
567 // * NASSX – Maximum number of assemblies across the core horizontally in full symmetry
568 // * NASSY – Maximum number of assemblies down the core vertically in full symmetry
569 // * NPIN – Maximum number of fuel pins across a fuel assembly in the core. Assemblies are
570 // assumed to be symmetric.
571 // * NAX – Number of axial levels edited in the fuel
572 // * NASS – Total number of fuel assemblies in the problem considering the symmetry of the
573 // calculation.
574 // --------------------------------------------------------------------------
575
576 this->ZCoordinates = this->ReadDataSet("/CORE", "axial_mesh");
577 this->NAX = this->Dimensions[0] - 1;
578 this->ZCoordinates->Delete(); // Let the smart pointer hold the only ref
579
580 this->CoreMap = this->ReadDataSet("/CORE", "core_map");
581 this->NASSX = this->Dimensions[0];
582 this->NASSY = this->Dimensions[1];
583 this->CoreMap->Delete(); // Let the smart pointer hold the only ref
584
585 dataSource = this->ReadDataSet("/CORE", "core_sym");
586 this->SYMMETRY = dataSource->GetTuple1(0);
587 dataSource->Delete(); // Just needed to extract the single value
588
589 // ------------------------------------------
590 // Extract pin information
591 // ------------------------------------------
592 std::vector<std::string> names;
593 this->AddDataSetNamesWithDimension("/CORE", 4, names);
594 for (const std::string& name : names)
595 {
596 dataSource = this->ReadDataSet("/CORE", name.c_str());
597 this->NPIN = this->Dimensions[0];
598 this->NASS = this->Dimensions[3];
599
600 outputCellArray = this->CreatePinFieldArray(dataSource);
601 outputCellArray->SetName(dataSource->GetName());
602 this->CoreCellData.push_back(outputCellArray);
603 outputCellArray->Delete();
604 dataSource->Delete();
605 }
606 // ------------------------------------------
607
608 // DEBUG: std::cout << "============================" << std::endl;
609 // DEBUG: std::cout << "NASSX " << this->NASSX << std::endl;
610 // DEBUG: std::cout << "NASSY " << this->NASSY << std::endl;
611 // DEBUG: std::cout << "NPIN " << this->NPIN << std::endl;
612 // DEBUG: std::cout << "NAX " << this->NAX << std::endl;
613 // DEBUG: std::cout << "NASS " << this->NASS << std::endl;
614 // DEBUG: std::cout << "SYMMETRY " << this->SYMMETRY << std::endl;
615 // DEBUG: std::cout << "============================" << std::endl;
616
617 // ------------------------------------------
618 // X/Y Coordinates
619 // ------------------------------------------
620
621 float axialStep = this->APITCH / (double)NPIN;
622 this->XCoordinates->SetNumberOfTuples(NASSX * NPIN + 1);
623 for (vtkIdType idx = 0; idx < this->XCoordinates->GetNumberOfTuples(); idx++)
624 {
625 this->XCoordinates->SetTuple1(idx, ((float)idx) * axialStep);
626 }
627
628 this->YCoordinates->SetNumberOfTuples(NASSY * NPIN + 1);
629 for (vtkIdType idx = 0; idx < this->YCoordinates->GetNumberOfTuples(); idx++)
630 {
631 this->YCoordinates->SetTuple1(idx, ((float)idx) * axialStep);
632 }
633
634 // ------------------------------------------
635 // Fill cellData from core information
636 // ------------------------------------------
637
638 outputCellArray = this->CoreMap->NewInstance();
639 outputCellArray->SetNumberOfTuples(
640 this->NASSX * this->NPIN * this->NASSY * this->NPIN * this->NAX);
641 outputCellArray->SetName("AssemblyID");
642
643 for (hsize_t sj = 0; sj < this->NASSY; sj++)
644 {
645 for (hsize_t si = 0; si < this->NASSX; si++)
646 {
647 vtkIdType srcIdx = si * this->NASSX + sj; // Fortran ordering
648 vtkIdType dstOffset = si * this->NPIN + sj * this->NASSX * this->NPIN * this->NPIN;
649 for (hsize_t dk = 0; dk < this->NAX; dk++)
650 {
651 for (hsize_t dj = 0; dj < this->NPIN; dj++)
652 {
653 for (hsize_t di = 0; di < this->NPIN; di++)
654 {
655 outputCellArray->SetTuple1(dstOffset + di + (dj * this->NASSX * this->NPIN) +
656 (dk * this->NASSX * this->NPIN * this->NASSY * this->NPIN),
657 this->CoreMap->GetTuple1(srcIdx));
658 }
659 }
660 }
661 }
662 }
663 this->CoreCellData.push_back(outputCellArray);
664 outputCellArray->Delete();
665 }
666
667 //----------------------------------------------------------------------------
668
InitializeWithCoreData(vtkRectilinearGrid * output)669 void InitializeWithCoreData(vtkRectilinearGrid* output)
670 {
671 // NoOp if already loaded
672 this->ReadCore();
673
674 output->SetDimensions(
675 this->NASSX * this->NPIN + 1, this->NASSY * this->NPIN + 1, this->NAX + 1);
676 output->SetXCoordinates(this->XCoordinates);
677 output->SetYCoordinates(this->YCoordinates);
678 output->SetZCoordinates(this->ZCoordinates);
679
680 for (const vtkSmartPointer<vtkDataArray>& cellArray : this->CoreCellData)
681 {
682 output->GetCellData()->AddArray(cellArray);
683 }
684 }
685
686 //----------------------------------------------------------------------------
687
AddStateData(vtkRectilinearGrid * output,vtkIdType timestep)688 void AddStateData(vtkRectilinearGrid* output, vtkIdType timestep)
689 {
690 if (this->FileId == -1)
691 {
692 return;
693 }
694
695 // --------------------------------------------------------------------------
696 // STATE_0001 Group
697 // STATE_0001/crit_boron Dataset {1}
698 // STATE_0001/detector_response Dataset {49, 18} <------ SKIP/FIXME
699 // STATE_0001/exposure Dataset {1}
700 // STATE_0001/keff Dataset {1}
701 // STATE_0001/pin_cladtemps Dataset {17, 17, 49, 56}
702 // STATE_0001/pin_fueltemps Dataset {17, 17, 49, 56}
703 // STATE_0001/pin_moddens Dataset {17, 17, 49, 56}
704 // STATE_0001/pin_modtemps Dataset {17, 17, 49, 56}
705 // STATE_0001/pin_powers Dataset {17, 17, 49, 56}
706 // --------------------------------------------------------------------------
707 std::ostringstream stateGroupName;
708 stateGroupName << "/STATE_" << std::setw(4) << std::setfill('0') << timestep;
709 vtkDataArray* dataSource;
710 vtkDataArray* outputCellArray;
711
712 // Open group
713 hid_t groupId = -1;
714 if ((groupId = H5Gopen(this->FileId, stateGroupName.str().c_str())) < 0)
715 {
716 vtkErrorWithObjectMacro(this->Owner, "Can't open Group " << stateGroupName.str().c_str());
717 return;
718 }
719
720 H5G_info_t groupInfo;
721 int status = H5Gget_info(groupId, &groupInfo);
722 if (status < 0)
723 {
724 vtkErrorWithObjectMacro(
725 this->Owner, "Can't get group info for " << stateGroupName.str().c_str());
726 return;
727 }
728
729 // Extract dataset names
730 std::vector<std::string> datasetNames;
731 char datasetName[DATASET_NAME_MAX_SIZE];
732 for (hsize_t idx = 0; idx < groupInfo.nlinks; idx++)
733 {
734 H5Lget_name_by_idx(groupId, ".", H5_INDEX_NAME, H5_ITER_INC, idx, datasetName,
735 DATASET_NAME_MAX_SIZE, H5P_DEFAULT);
736 datasetNames.push_back(datasetName);
737 }
738 H5Gclose(groupId);
739
740 // Start processing datasets
741 for (const std::string& dsName : datasetNames)
742 {
743 if (!(this->CellDataArraySelection->ArrayExists(dsName.c_str()) ||
744 this->FieldDataArraySelection->ArrayExists(dsName.c_str())) ||
745 !(this->CellDataArraySelection->ArrayIsEnabled(dsName.c_str()) ||
746 this->FieldDataArraySelection->ArrayIsEnabled(dsName.c_str())))
747 {
748 continue;
749 }
750 dataSource = this->ReadDataSet(stateGroupName.str().c_str(), dsName.c_str());
751 if (dataSource)
752 {
753 if (this->NumberOfDimensions == 4 && this->Dimensions[0] == this->NPIN &&
754 this->Dimensions[1] == this->NPIN && this->Dimensions[2] == this->NAX &&
755 this->Dimensions[3] == this->NASS)
756 {
757 outputCellArray = this->CreatePinFieldArray(dataSource);
758 outputCellArray->SetName(dsName.c_str());
759 output->GetCellData()->AddArray(outputCellArray);
760 outputCellArray->Delete();
761 }
762 else if (this->NumberOfDimensions == 1 && this->Dimensions[0] == 1)
763 {
764 // Add fields add FieldData
765 output->GetFieldData()->AddArray(dataSource);
766 }
767 else
768 {
769 std::ostringstream message;
770 message << "Invalid dimensions: ";
771 for (hsize_t i = 0; i < this->NumberOfDimensions; i++)
772 {
773 message << this->Dimensions[i] << " ";
774 }
775 vtkDebugWithObjectMacro(this->Owner, << message.str());
776 }
777 dataSource->Delete();
778 }
779 }
780 }
781
782 //----------------------------------------------------------------------------
783
784 vtkNew<vtkDataArraySelection> CellDataArraySelection;
785 vtkNew<vtkDataArraySelection> FieldDataArraySelection;
786
787 private:
788 hid_t FileId;
789 std::string FileName;
790
791 hsize_t NumberOfDimensions;
792 hsize_t Dimensions[VERA_MAX_DIMENSION];
793
794 bool NeedCoreProcessing;
795
796 double APITCH;
797 hsize_t NASSX;
798 hsize_t NASSY;
799 hsize_t NAX;
800 hsize_t NPIN;
801 hsize_t NASS;
802 vtkIdType SYMMETRY;
803 vtkIdType NUMBER_OF_STATES;
804
805 vtkNew<vtkFloatArray> XCoordinates;
806 vtkNew<vtkFloatArray> YCoordinates;
807
808 vtkObject* Owner;
809 vtkSmartPointer<vtkDataArray> ZCoordinates;
810 vtkSmartPointer<vtkDataArray> CoreMap;
811
812 std::vector<vtkSmartPointer<vtkDataArray> > CoreCellData;
813 };
814 //*****************************************************************************
815
816 vtkStandardNewMacro(vtkVeraOutReader);
817
818 //----------------------------------------------------------------------------
819 // Constructor for vtkVeraOutReader
820 //----------------------------------------------------------------------------
vtkVeraOutReader()821 vtkVeraOutReader::vtkVeraOutReader()
822 {
823 this->FileName = nullptr;
824 this->NumberOfTimeSteps = 0;
825 this->TimeSteps.clear();
826 this->SetNumberOfInputPorts(0);
827 this->SetNumberOfOutputPorts(1);
828 this->Internal = new Internals(this);
829 }
830
831 //----------------------------------------------------------------------------
832 // Destructor for vtkVeraOutReader
833 //----------------------------------------------------------------------------
~vtkVeraOutReader()834 vtkVeraOutReader::~vtkVeraOutReader()
835 {
836 this->SetFileName(nullptr);
837 delete this->Internal;
838 this->Internal = nullptr;
839 }
840
841 //----------------------------------------------------------------------------
842 // Verify that the file exists, get dimension sizes and variables
843 //----------------------------------------------------------------------------
RequestInformation(vtkInformation * reqInfo,vtkInformationVector ** inVector,vtkInformationVector * outVector)844 int vtkVeraOutReader::RequestInformation(
845 vtkInformation* reqInfo, vtkInformationVector** inVector, vtkInformationVector* outVector)
846 {
847 vtkDebugMacro(<< "In vtkVeraOutReader::RequestInformation" << endl);
848 if (!this->Superclass::RequestInformation(reqInfo, inVector, outVector))
849 return 0;
850
851 if (!this->FileName || this->FileName[0] == '\0')
852 {
853 vtkErrorMacro("No filename specified");
854 return 0;
855 }
856
857 vtkDebugMacro(<< "In vtkVeraOutReader::RequestInformation read filename okay" << endl);
858 vtkInformation* outInfo = outVector->GetInformationObject(0);
859
860 this->NumberOfTimeSteps = 0;
861 this->Internal->SetFileName(this->FileName);
862 if (this->Internal->OpenFile())
863 {
864 this->NumberOfTimeSteps = this->Internal->GetNumberOfTimeSteps();
865 this->Internal->LoadMetaData();
866 this->Internal->CloseFile();
867 }
868
869 this->TimeSteps.resize(this->NumberOfTimeSteps);
870 std::iota(this->TimeSteps.begin(), this->TimeSteps.end(), 1.0);
871
872 outInfo->Set(
873 vtkStreamingDemandDrivenPipeline::TIME_STEPS(), &this->TimeSteps[0], this->NumberOfTimeSteps);
874
875 double tRange[2];
876 tRange[0] = this->NumberOfTimeSteps ? this->TimeSteps[0] : 0;
877 tRange[1] = this->NumberOfTimeSteps ? this->TimeSteps[this->NumberOfTimeSteps - 1] : 0;
878 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), tRange, 2);
879
880 return 1;
881 }
882
883 //----------------------------------------------------------------------------
884 // Default method: Data is read into a vtkUnstructuredGrid
885 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (reqInfo),vtkInformationVector ** vtkNotUsed (inVector),vtkInformationVector * outVector)886 int vtkVeraOutReader::RequestData(vtkInformation* vtkNotUsed(reqInfo),
887 vtkInformationVector** vtkNotUsed(inVector), vtkInformationVector* outVector)
888 {
889 if (!this->FileName || this->FileName[0] == '\0')
890 {
891 vtkErrorMacro("No filename specified");
892 return 0;
893 }
894
895 vtkDebugMacro(<< "In vtkVeraOutReader::RequestData" << endl);
896 vtkInformation* outInfo = outVector->GetInformationObject(0);
897 vtkRectilinearGrid* output =
898 vtkRectilinearGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
899
900 // --------------------------------------------------------------------------
901 // Time / State handling
902 // --------------------------------------------------------------------------
903
904 vtkIdType requestedTimeStep = 0;
905 auto timeKey = vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP();
906
907 if (outInfo->Has(timeKey))
908 {
909 requestedTimeStep = outInfo->Get(timeKey);
910 }
911
912 // --------------------------------------------------------------------------
913 // Data handling
914 // --------------------------------------------------------------------------
915
916 this->Internal->SetFileName(this->FileName);
917 if (this->Internal->OpenFile())
918 {
919 this->Internal->InitializeWithCoreData(output);
920 this->Internal->AddStateData(output, requestedTimeStep);
921 this->Internal->CloseFile();
922 }
923
924 // --------------------------------------------------------------------------
925
926 vtkDebugMacro(<< "Out vtkVeraOutReader::RequestData" << endl);
927
928 return 1;
929 }
930
931 //----------------------------------------------------------------------------
932 // Print self.
933 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)934 void vtkVeraOutReader::PrintSelf(ostream& os, vtkIndent indent)
935 {
936 this->Superclass::PrintSelf(os, indent);
937 os << indent << "FileName: " << (this->FileName ? this->FileName : "NULL") << "\n";
938 }
939
940 //----------------------------------------------------------------------------
941 // Cell Array Selection
942 //----------------------------------------------------------------------------
943
GetCellDataArraySelection() const944 vtkDataArraySelection* vtkVeraOutReader::GetCellDataArraySelection() const
945 {
946 return this->Internal->CellDataArraySelection;
947 }
948
949 //----------------------------------------------------------------------------
950 // Field Array Selection
951 //----------------------------------------------------------------------------
952
GetFieldDataArraySelection() const953 vtkDataArraySelection* vtkVeraOutReader::GetFieldDataArraySelection() const
954 {
955 return this->Internal->FieldDataArraySelection;
956 }
957
958 //----------------------------------------------------------------------------
959 // MTime
960 //----------------------------------------------------------------------------
961
GetMTime()962 vtkMTimeType vtkVeraOutReader::GetMTime()
963 {
964 vtkMTimeType mTime = this->Superclass::GetMTime();
965 vtkMTimeType cellMTime = this->Internal->CellDataArraySelection->GetMTime();
966 vtkMTimeType fieldMTime = this->Internal->FieldDataArraySelection->GetMTime();
967
968 mTime = ( cellMTime > mTime ? cellMTime : mTime );
969 mTime = ( fieldMTime > mTime ? fieldMTime : mTime );
970
971 return mTime;
972 }
973