1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkAMReXGridReaderInternal.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 #include "vtkAMReXGridReaderInternal.h"
17 #include "vtkByteSwap.h"
18 #include "vtkCellData.h"
19 #include "vtkDoubleArray.h"
20 #include "vtkFloatArray.h"
21 #include "vtkIndent.h"
22 #include "vtkObject.h"
23 #include "vtkSetGet.h"
24
25 #include <sstream>
26 #include <vector>
27 #include <vtksys/FStream.hxx>
28
29 namespace
30 {
ReadFile(const std::string & filename)31 std::string ReadFile(const std::string& filename)
32 {
33 std::string contents;
34 vtksys::ifstream stream(filename.c_str(), std::ios::binary);
35 if (stream)
36 {
37 stream.seekg(0, std::ios::end);
38 int flength = static_cast<int>(stream.tellg());
39 stream.seekg(0, std::ios::beg);
40 std::vector<char> data(flength + 1 + (flength + 1) % 8); // padded for better alignment.
41 stream.read(data.data(), flength);
42 data[flength] = '\0';
43 contents = data.data();
44 }
45 return contents;
46 }
47 }
48
49 RealDescriptor::RealDescriptor() = default;
50
RealDescriptor(const long * fr_,const int * ord_,int ordl_)51 RealDescriptor::RealDescriptor(const long* fr_, const int* ord_, int ordl_)
52 : fr(fr_, fr_ + 8)
53 , ord(ord_, ord_ + ordl_)
54 {
55 }
56
format() const57 const long* RealDescriptor::format() const&
58 {
59 return &fr[0];
60 }
61
formatarray() const62 const std::vector<long>& RealDescriptor::formatarray() const&
63 {
64 return fr;
65 }
66
order() const67 const int* RealDescriptor::order() const&
68 {
69 return &ord[0];
70 }
71
orderarray() const72 const std::vector<int>& RealDescriptor::orderarray() const&
73 {
74 return ord;
75 }
76
numBytes() const77 int RealDescriptor::numBytes() const
78 {
79 return (fr[0] + 7) >> 3;
80 }
81
operator ==(const RealDescriptor & rd) const82 bool RealDescriptor::operator==(const RealDescriptor& rd) const
83 {
84 return fr == rd.fr && ord == rd.ord;
85 }
86
87 #define AMREX_PRINT(os, indent, var) os << indent << #var << ": " << var << endl
88
vtkAMReXGridHeader()89 vtkAMReXGridHeader::vtkAMReXGridHeader()
90 : variableNamesSize(0)
91 , dim(0)
92 , time(0)
93 , finestLevel(0)
94 , geometryCoord(0)
95 , magicZero(0)
96 , debugHeader(false)
97 {
98 }
99
PrintSelf(std::ostream & os,vtkIndent indent)100 void vtkAMReXGridHeader::PrintSelf(std::ostream& os, vtkIndent indent)
101 {
102 this->PrintSelfGenericHeader(os, indent);
103 }
104
PrintSelfGenericHeader(ostream & os,vtkIndent indent)105 void vtkAMReXGridHeader::PrintSelfGenericHeader(ostream& os, vtkIndent indent)
106 {
107 AMREX_PRINT(os, indent, versionName);
108 AMREX_PRINT(os, indent, variableNamesSize);
109 os << indent << "variableNames: " << endl;
110 for (const auto& name : this->variableNames)
111 {
112 os << indent.GetNextIndent() << name << endl;
113 }
114 AMREX_PRINT(os, indent, dim);
115 AMREX_PRINT(os, indent, time);
116 AMREX_PRINT(os, indent, finestLevel);
117 os << indent << "problemDomainLoEnd: " << endl << indent.GetNextIndent();
118 for (const auto& lo : this->problemDomainLoEnd)
119 {
120 os << lo << " ";
121 }
122 os << endl;
123 os << indent << "problemDomainHiEnd: " << endl << indent.GetNextIndent();
124 for (const auto& hi : this->problemDomainHiEnd)
125 {
126 os << hi << " ";
127 }
128 os << endl;
129 os << indent << "refinementRatio: " << endl << indent.GetNextIndent();
130 for (const auto& ratio : this->refinementRatio)
131 {
132 os << ratio << " ";
133 }
134 os << endl;
135 os << indent << "levelDomains: " << endl << indent.GetNextIndent();
136 for (int level = 0; level < this->finestLevel + 1; ++level)
137 {
138 os << "(";
139 os << "(";
140 for (int space = 0; space < this->dim; ++space)
141 {
142 os << this->levelDomains[level][0][space];
143 if (space < (this->dim - 1))
144 os << ",";
145 }
146 os << ") ";
147 os << "(";
148 for (int space = 0; space < this->dim; ++space)
149 {
150 os << this->levelDomains[level][1][space];
151 if (space < (this->dim - 1))
152 os << ",";
153 }
154 os << ") ";
155 os << "(";
156 for (int space = 0; space < this->dim; ++space)
157 {
158 os << this->levelDomains[level][2][space];
159 if (space < (this->dim - 1))
160 os << ",";
161 }
162 os << ")";
163 if (level < this->finestLevel)
164 os << ") ";
165 else
166 os << ")";
167 }
168 os << endl;
169 os << indent << "levelSteps: " << endl << indent.GetNextIndent();
170 for (const auto& steps : this->levelSteps)
171 {
172 os << steps << " ";
173 }
174 os << endl;
175 os << indent << "cellSize: " << endl << indent.GetNextIndent();
176 for (int level = 0; level < this->finestLevel + 1; ++level)
177 {
178 for (int space = 0; space < this->dim; ++space)
179 {
180 os << this->cellSize[level][space];
181 }
182 if (level < this->finestLevel)
183 os << endl << indent.GetNextIndent();
184 else
185 os << endl;
186 }
187 AMREX_PRINT(os, indent, geometryCoord);
188 AMREX_PRINT(os, indent, magicZero);
189 os << indent << "levelCells: " << endl << indent.GetNextIndent();
190 for (int level = 0; level <= this->finestLevel; ++level)
191 {
192 os << level << " " << this->levelSize[level] << " " << this->time << endl
193 << indent.GetNextIndent();
194 os << this->levelSteps[level] << endl << indent.GetNextIndent();
195 for (int size = 0; size < this->levelSize[level]; ++size)
196 {
197 for (int space = 0; space < this->dim; ++space)
198 {
199 for (int bound = 0; bound <= 1; ++bound)
200 {
201 os << this->levelCells[level][size][space][bound] << " ";
202 }
203 os << endl << indent.GetNextIndent();
204 }
205 }
206 os << this->levelPrefix[level] << "/" << this->multiFabPrefix[level] << endl
207 << indent.GetNextIndent();
208 }
209 os << "Generic Header Complete" << endl;
210 }
211
Parse(const std::string & headerData)212 bool vtkAMReXGridHeader::Parse(const std::string& headerData)
213 {
214 this->ParseGenericHeader(headerData);
215 if (debugHeader)
216 {
217 std::ostream* os;
218 os = &std::cout;
219 vtkIndent indent;
220 this->PrintSelf(*os, indent);
221 }
222 return true;
223 }
224
SetVectorNamePrefix(const std::string & prefix)225 void vtkAMReXGridHeader::SetVectorNamePrefix(const std::string& prefix)
226 {
227 this->vectorNamePrefix = prefix;
228 }
229
SetNameDelimiter(const char delim)230 void vtkAMReXGridHeader::SetNameDelimiter(const char delim)
231 {
232 this->nameDelim = delim;
233 }
234
GetBaseVariableName(const std::string & name)235 std::string vtkAMReXGridHeader::GetBaseVariableName(const std::string& name)
236 {
237 std::string baseName = name;
238 const std::size_t pos = baseName.find_first_of(this->nameDelim);
239 const std::string prefix = baseName.substr(0, pos);
240 if (prefix == this->vectorNamePrefix)
241 {
242 baseName = baseName.substr(pos + 1);
243 std::size_t posSuffix = baseName.find_last_of(this->nameDelim);
244 baseName = baseName.substr(0, posSuffix);
245 }
246 return baseName;
247 }
248
CheckComponent(const std::string & name)249 int vtkAMReXGridHeader::CheckComponent(const std::string& name)
250 {
251 const std::size_t pos = name.find_last_of(this->nameDelim);
252 // we expect to use the character just past pos
253 // so we don't want to accidentially jump outside the string
254 if (pos > name.size() - 1)
255 {
256 return -1;
257 }
258 const std::string suffix = name.substr(pos + 1);
259 if (suffix == "x")
260 {
261 return 0;
262 }
263 if (suffix == "y")
264 {
265 return 1;
266 }
267 if (suffix == "z")
268 {
269 return 2;
270 }
271 return -1;
272 }
273
ParseGenericHeader(const std::string & headerData)274 bool vtkAMReXGridHeader::ParseGenericHeader(const std::string& headerData)
275 {
276 char c;
277 std::istringstream hstream(headerData);
278 hstream >> this->versionName;
279 if (this->versionName.empty())
280 {
281 vtkGenericWarningMacro("Failed to read versionName string.");
282 return false;
283 }
284 hstream >> this->variableNamesSize;
285 this->variableNames.resize(this->variableNamesSize);
286 for (int cc = 0; cc < this->variableNamesSize; ++cc)
287 {
288 hstream >> this->variableNames[cc];
289
290 // build vector name map
291 std::string baseName = this->GetBaseVariableName(this->variableNames[cc]);
292 int component = this->CheckComponent(this->variableNames[cc]);
293 // resize variable map for scalar or vector
294 parsedVariableNames[baseName].resize(component < 0 ? 1 : 3);
295 component = (component < 0) ? 0 : component;
296 parsedVariableNames[baseName][component] = cc;
297 }
298 hstream >> this->dim;
299 if (this->dim != 1 && this->dim != 2 && this->dim != 3)
300 {
301 vtkGenericWarningMacro("dim must be 1, 2, or 3.");
302 return false;
303 }
304 hstream >> this->time;
305 hstream >> this->finestLevel;
306 if (this->finestLevel < 0)
307 {
308 vtkGenericWarningMacro("finestLevel must be >= 0");
309 return false;
310 }
311 this->problemDomainLoEnd.resize(this->dim);
312 for (int cc = 0; cc < this->dim; ++cc)
313 {
314 hstream >> this->problemDomainLoEnd[cc];
315 }
316 this->problemDomainHiEnd.resize(this->dim);
317 for (int cc = 0; cc < this->dim; ++cc)
318 {
319 hstream >> this->problemDomainHiEnd[cc];
320 }
321 this->refinementRatio.resize(this->finestLevel);
322 for (int cc = 0; cc < this->finestLevel; ++cc)
323 {
324 hstream >> this->refinementRatio[cc];
325 }
326 this->levelDomains.resize(this->finestLevel + 1);
327 for (int cc = 0; cc <= this->finestLevel; ++cc)
328 {
329 this->levelDomains[cc].resize(3);
330 for (int dd = 0; dd < 3; ++dd)
331 {
332 this->levelDomains[cc][dd].resize(this->dim);
333 }
334 }
335 for (int cc = 0; cc <= this->finestLevel; ++cc)
336 {
337 hstream >> c; // read '('
338 for (int dd = 0; dd < 3; ++dd)
339 {
340 hstream >> c; // read '('
341 for (int ee = 0; ee < this->dim; ++ee)
342 {
343 hstream >> this->levelDomains[cc][dd][ee];
344 if (ee < (this->dim - 1))
345 {
346 hstream >> c; // read ','
347 }
348 else
349 {
350 hstream >> c; // read ')'
351 }
352 }
353 }
354 hstream >> c; // read ')'
355 }
356 this->levelSteps.resize(this->finestLevel + 1);
357 for (int cc = 0; cc < this->finestLevel + 1; ++cc)
358 {
359 hstream >> this->levelSteps[cc];
360 }
361 this->cellSize.resize(this->finestLevel + 1);
362 for (int cc = 0; cc <= this->finestLevel; ++cc)
363 {
364 this->cellSize[cc].resize(this->dim);
365 }
366 for (int cc = 0; cc <= this->finestLevel; ++cc)
367 {
368 for (int dd = 0; dd < this->dim; ++dd)
369 {
370 hstream >> this->cellSize[cc][dd];
371 }
372 }
373 hstream >> this->geometryCoord;
374 hstream >> this->magicZero;
375 int tmpLevel;
376 this->levelSize.resize(this->finestLevel + 1);
377 double tmpTime;
378 int tmpLevelSteps;
379 this->levelCells.resize(this->finestLevel + 1);
380 std::string pathName;
381 std::string deliminator = "/";
382 this->levelPrefix.resize(this->finestLevel + 1);
383 this->multiFabPrefix.resize(this->finestLevel + 1);
384 for (int level = 0; level <= this->finestLevel; ++level)
385 {
386 hstream >> tmpLevel;
387 hstream >> this->levelSize[level];
388 hstream >> tmpTime;
389 hstream >> tmpLevelSteps;
390 this->levelCells[level].resize(this->levelSize[level]);
391 for (int size = 0; size < this->levelSize[level]; ++size)
392 {
393 this->levelCells[level][size].resize(this->dim);
394 for (int space = 0; space < this->dim; ++space)
395 {
396 this->levelCells[level][size][space].resize(2);
397 }
398 }
399 for (int size = 0; size < this->levelSize[level]; ++size)
400 {
401 for (int space = 0; space < this->dim; ++space)
402 {
403 for (int bound = 0; bound <= 1; ++bound)
404 {
405 hstream >> this->levelCells[level][size][space][bound];
406 }
407 }
408 }
409 hstream >> pathName;
410 size_t pos = 0;
411 pos = pathName.find(deliminator);
412 this->levelPrefix[level] = pathName.substr(0, pos);
413 pathName.erase(0, pos + deliminator.length());
414 this->multiFabPrefix[level] = pathName.substr(0, pathName.length());
415 }
416 return true;
417 }
418
vtkAMReXGridLevelHeader()419 vtkAMReXGridLevelHeader::vtkAMReXGridLevelHeader()
420 : level(0)
421 , dim(0)
422 , levelVersion(0)
423 , levelHow(0)
424 , levelNumberOfComponents(0)
425 , levelNumberOfGhostCells(0)
426 , levelBoxArraySize(0)
427 , levelMagicZero(0)
428 , levelNumberOfFABOnDisk(0)
429 , levelRealNumberOfBytes(0)
430 , levelRealOrder(0)
431 , debugLevelHeader(false)
432 {
433 }
434
PrintSelf(ostream & os,vtkIndent indent)435 void vtkAMReXGridLevelHeader::PrintSelf(ostream& os, vtkIndent indent)
436 {
437 this->PrintSelfLevelHeader(os, indent);
438 }
439
PrintSelfLevelHeader(ostream & os,vtkIndent indent)440 void vtkAMReXGridLevelHeader::PrintSelfLevelHeader(ostream& os, vtkIndent indent)
441 {
442 AMREX_PRINT(os, indent, level);
443 AMREX_PRINT(os, indent, levelVersion);
444 AMREX_PRINT(os, indent, levelHow);
445 AMREX_PRINT(os, indent, levelNumberOfComponents);
446 AMREX_PRINT(os, indent, levelNumberOfGhostCells);
447 os << "BoxArray Size and MagicZero:" << endl << indent.GetNextIndent();
448 os << "(" << this->levelBoxArraySize << " " << this->levelMagicZero << endl
449 << indent.GetNextIndent();
450 os << indent << "levelDomains: " << endl << indent.GetNextIndent();
451 for (int cc = 0; cc < this->levelBoxArraySize; ++cc)
452 {
453 os << "(";
454 os << "(";
455 for (int space = 0; space < this->dim; ++space)
456 {
457 os << this->levelBoxArrays[cc][0][space];
458 if (space < (this->dim - 1))
459 os << ",";
460 }
461 os << ") ";
462 os << "(";
463 for (int space = 0; space < this->dim; ++space)
464 {
465 os << this->levelBoxArrays[cc][1][space];
466 if (space < (this->dim - 1))
467 os << ",";
468 }
469 os << ") ";
470 os << "(";
471 for (int space = 0; space < this->dim; ++space)
472 {
473 os << this->levelBoxArrays[cc][2][space];
474 if (space < (this->dim - 1))
475 os << ",";
476 }
477 os << ")";
478 if (cc < (this->levelBoxArraySize - 1))
479 os << ")" << endl << indent.GetNextIndent();
480 else
481 os << ")" << endl;
482 }
483 os << ")" << endl;
484 AMREX_PRINT(os, indent, levelNumberOfFABOnDisk);
485 os << indent << "FABsOnDisk: " << endl << indent.GetNextIndent();
486 for (int cc = 0; cc < this->levelNumberOfFABOnDisk; ++cc)
487 {
488 os << this->levelFabOnDiskPrefix << ' ' << this->levelFABFile[cc] << ' '
489 << this->levelFileOffset[cc] << endl
490 << indent.GetNextIndent();
491 }
492 os << endl;
493 if (this->levelVersion == Version_v1 || this->levelVersion == NoFabHeaderMinMax_v1)
494 {
495 std::ios::fmtflags oflags = os.flags();
496 os.setf(std::ios::floatfield, std::ios::scientific);
497 int oldPrec(os.precision(16));
498 os << indent << "Minimums and Maximums of FABs: " << endl << indent.GetNextIndent();
499 os << this->levelNumberOfFABOnDisk << "," << this->levelNumberOfComponents << endl
500 << indent.GetNextIndent();
501 for (int cc = 0; cc < this->levelNumberOfFABOnDisk; ++cc)
502 {
503 for (int dd = 0; dd < this->levelNumberOfComponents; ++dd)
504 {
505 os << this->levelMinimumsFAB[cc][dd] << "," << endl << indent.GetNextIndent();
506 }
507 }
508 os << endl << indent.GetNextIndent();
509 os << this->levelNumberOfFABOnDisk << "," << this->levelNumberOfComponents << endl
510 << indent.GetNextIndent();
511 for (int cc = 0; cc < this->levelNumberOfFABOnDisk; ++cc)
512 {
513 for (int dd = 0; dd < this->levelNumberOfComponents; ++dd)
514 {
515 os << this->levelMaximumsFAB[cc][dd];
516 if (cc < (this->levelNumberOfFABOnDisk - 1))
517 os << "," << endl << indent.GetNextIndent();
518 else
519 os << endl << indent.GetNextIndent();
520 }
521 }
522 os << endl << indent.GetNextIndent();
523 os.flags(oflags);
524 os.precision(oldPrec);
525 }
526 if (this->levelVersion == NoFabHeaderFAMinMax_v1)
527 {
528 os << indent << "Minimums and Maximums of FABArray: " << endl << indent.GetNextIndent();
529 for (int cc = 0; cc < this->levelNumberOfComponents; ++cc)
530 {
531 os << this->levelFABArrayMinimum[cc] << ",";
532 }
533 os << endl << indent.GetNextIndent();
534 for (int cc = 0; cc < this->levelNumberOfComponents; ++cc)
535 {
536 os << this->levelFABArrayMaximum[cc] << ",";
537 }
538 os << endl << indent.GetNextIndent();
539 }
540 if (this->levelVersion == NoFabHeader_v1 || this->levelVersion == NoFabHeaderMinMax_v1 ||
541 this->levelVersion == NoFabHeaderFAMinMax_v1)
542 {
543 os << indent << "Real Format: " << endl << indent.GetNextIndent();
544 os << "(" << this->levelRealNumberOfBytes << "," << this->levelRealOrder << ")" << endl
545 << indent.GetNextIndent();
546 }
547 os << "Level " << this->level << " Header Complete" << endl;
548 }
549
Parse(int _level,int _dim,const std::string & headerData)550 bool vtkAMReXGridLevelHeader::Parse(int _level, int _dim, const std::string& headerData)
551 {
552 this->ParseLevelHeader(_level, _dim, headerData);
553 if (debugLevelHeader)
554 {
555 std::ostream* os;
556 os = &std::cout;
557 vtkIndent indent;
558 this->PrintSelf(*os, indent);
559 }
560 return true;
561 }
562
ParseLevelHeader(int _level,int _dim,const std::string & headerData)563 bool vtkAMReXGridLevelHeader::ParseLevelHeader(int _level, int _dim, const std::string& headerData)
564 {
565 std::istringstream hstream(headerData);
566 this->level = _level;
567 this->dim = _dim;
568 hstream >> this->levelVersion;
569 hstream >> this->levelHow;
570 hstream >> this->levelNumberOfComponents;
571 hstream >> this->levelNumberOfGhostCells;
572 char c;
573 hstream >> c; // read '(' Begin BoxArray writeon()
574 hstream >> this->levelBoxArraySize;
575 hstream >> this->levelMagicZero;
576 this->levelBoxArrays.resize(this->levelBoxArraySize);
577 for (int cc = 0; cc < this->levelBoxArraySize; ++cc)
578 {
579 this->levelBoxArrays[cc].resize(3);
580 for (int dd = 0; dd < 3; ++dd)
581 {
582 this->levelBoxArrays[cc][dd].resize(this->dim);
583 }
584 }
585 for (int cc = 0; cc < this->levelBoxArraySize; ++cc)
586 {
587 hstream >> c; // read '('
588 for (int dd = 0; dd < 3; ++dd)
589 {
590 hstream >> c; // read '('
591 for (int ee = 0; ee < this->dim; ++ee)
592 {
593 hstream >> this->levelBoxArrays[cc][dd][ee];
594 if (ee < (this->dim - 1))
595 {
596 hstream >> c; // read ','
597 }
598 else
599 {
600 hstream >> c; // read ')'
601 }
602 }
603 }
604 hstream >> c; // read ')'
605 }
606 hstream >> c; // read ')' End BoxArray writeon()
607 hstream >> this->levelNumberOfFABOnDisk;
608 this->levelFABFile.resize(this->levelNumberOfFABOnDisk);
609 this->levelFileOffset.resize(this->levelNumberOfFABOnDisk);
610 for (int cc = 0; cc < this->levelNumberOfFABOnDisk; ++cc)
611 {
612 hstream >> this->levelFabOnDiskPrefix; // Prefix
613 hstream >> this->levelFABFile[cc]; // File
614 hstream >> this->levelFileOffset[cc]; // Offset
615 }
616 if (this->levelVersion == Version_v1 || this->levelVersion == NoFabHeaderMinMax_v1)
617 {
618 int tmpLevelNumberOfFABOnDisk;
619 int tmpLevelNumberOfComponents;
620 hstream >> tmpLevelNumberOfFABOnDisk;
621 hstream >> c; // read ','
622 hstream >> tmpLevelNumberOfComponents;
623 this->levelMinimumsFAB.resize(this->levelNumberOfFABOnDisk);
624 for (int cc = 0; cc < this->levelNumberOfFABOnDisk; ++cc)
625 {
626 this->levelMinimumsFAB[cc].resize(this->levelNumberOfComponents);
627 for (int dd = 0; dd < this->levelNumberOfComponents; ++dd)
628 {
629 hstream >> this->levelMinimumsFAB[cc][dd];
630 hstream >> c;
631 }
632 }
633 hstream >> tmpLevelNumberOfFABOnDisk;
634 hstream >> c; // read ','
635 hstream >> tmpLevelNumberOfComponents;
636 this->levelMaximumsFAB.resize(this->levelNumberOfFABOnDisk);
637 for (int cc = 0; cc < this->levelNumberOfFABOnDisk; ++cc)
638 {
639 this->levelMaximumsFAB[cc].resize(this->levelNumberOfComponents);
640 for (int dd = 0; dd < this->levelNumberOfComponents; ++dd)
641 {
642 hstream >> this->levelMaximumsFAB[cc][dd];
643 hstream >> c; // read ','
644 }
645 }
646 }
647 if (this->levelVersion == NoFabHeaderFAMinMax_v1)
648 {
649 this->levelFABArrayMinimum.resize(this->levelNumberOfComponents);
650 for (int cc = 0; cc < this->levelNumberOfComponents; ++cc)
651 {
652 hstream >> this->levelFABArrayMinimum[cc];
653 hstream >> c; // read ','
654 }
655 this->levelFABArrayMaximum.resize(this->levelNumberOfComponents);
656 for (int cc = 0; cc < this->levelNumberOfComponents; ++cc)
657 {
658 hstream >> this->levelFABArrayMaximum[cc];
659 hstream >> c; // read ','
660 }
661 }
662 if (this->levelVersion == NoFabHeader_v1 || this->levelVersion == NoFabHeaderMinMax_v1 ||
663 this->levelVersion == NoFabHeaderFAMinMax_v1)
664 {
665 hstream >> c; // read '('
666 hstream >> this->levelRealNumberOfBytes;
667 hstream >> c; // read ','
668 hstream >> this->levelRealOrder;
669 hstream >> c; // read ')'
670 }
671 return true;
672 }
673
vtkAMReXGridReaderInternal()674 vtkAMReXGridReaderInternal::vtkAMReXGridReaderInternal()
675 {
676 this->headersAreRead = false;
677 this->debugReader = false;
678 this->Header = nullptr;
679 }
680
~vtkAMReXGridReaderInternal()681 vtkAMReXGridReaderInternal::~vtkAMReXGridReaderInternal()
682 {
683 this->DestroyHeader();
684 this->DestroyLevelHeader();
685 }
686
DestroyHeader()687 void vtkAMReXGridReaderInternal::DestroyHeader()
688 {
689 delete this->Header;
690 this->Header = nullptr;
691 }
692
DestroyLevelHeader()693 void vtkAMReXGridReaderInternal::DestroyLevelHeader()
694 {
695 for (unsigned int i = 0; i < this->LevelHeader.size(); ++i)
696 {
697 delete this->LevelHeader[i];
698 this->LevelHeader[i] = nullptr;
699 }
700 }
701
PrintSelf(ostream & os,vtkIndent indent)702 void vtkAMReXGridReaderInternal::PrintSelf(ostream& os, vtkIndent indent)
703 {
704 os << indent << "FileName: " << this->FileName << endl;
705 if (this->Header)
706 {
707 os << indent << "Header: " << endl;
708 this->Header->PrintSelf(os, indent.GetNextIndent());
709 os << indent << "LevelHeader(s): " << endl;
710 for (int cc = 0; cc < this->Header->finestLevel + 1; ++cc)
711 {
712 this->LevelHeader[cc]->PrintSelfLevelHeader(os, indent.GetNextIndent());
713 }
714 }
715 else
716 {
717 os << indent << "Header: nullptr" << endl;
718 }
719 }
720
SetFileName(char * fName)721 void vtkAMReXGridReaderInternal::SetFileName(char* fName)
722 {
723 const std::string filename(fName == nullptr ? "" : fName);
724 this->FileName = filename;
725 this->headersAreRead = false;
726 }
727
728 //------------------------------------------------------------------------------
ReadMetaData()729 void vtkAMReXGridReaderInternal::ReadMetaData()
730 {
731 if (!this->headersAreRead)
732 {
733 if (!this->FileName.empty())
734 {
735 if (this->ReadHeader())
736 {
737 this->headersAreRead = this->ReadLevelHeader();
738 }
739 }
740 }
741 }
742
743 //------------------------------------------------------------------------------
ReadHeader()744 bool vtkAMReXGridReaderInternal::ReadHeader()
745 {
746 this->DestroyHeader();
747
748 const std::string headerFileName = this->FileName + "/Header";
749 const auto headerData = ::ReadFile(headerFileName);
750 if (headerData.empty())
751 {
752 return false;
753 }
754
755 auto headerPtr = new vtkAMReXGridHeader();
756 if (!headerPtr->Parse(headerData))
757 {
758 delete headerPtr;
759 return false;
760 }
761
762 this->Header = headerPtr;
763 return true;
764 }
765
766 //------------------------------------------------------------------------------
ReadLevelHeader()767 bool vtkAMReXGridReaderInternal::ReadLevelHeader()
768 {
769 this->DestroyLevelHeader();
770
771 this->LevelHeader.resize(this->Header->finestLevel + 1);
772
773 for (int level = 0; level <= this->Header->finestLevel; ++level)
774 {
775 const std::string levelHeaderFileName = this->FileName + "/" +
776 this->Header->levelPrefix[level] + "/" + this->Header->multiFabPrefix[level] + "_H";
777
778 const auto headerData = ::ReadFile(levelHeaderFileName);
779 if (headerData.empty())
780 {
781 return false;
782 }
783
784 auto headerPtr = new vtkAMReXGridLevelHeader();
785 if (!headerPtr->Parse(level, this->Header->dim, headerData))
786 {
787 delete headerPtr;
788 return false;
789 }
790
791 this->LevelHeader[level] = headerPtr;
792 }
793 return true;
794 }
795
GetNumberOfLevels()796 int vtkAMReXGridReaderInternal::GetNumberOfLevels()
797 {
798 return this->headersAreRead ? this->Header->finestLevel : -1;
799 }
800
GetBlockLevel(const int blockIdx)801 int vtkAMReXGridReaderInternal::GetBlockLevel(const int blockIdx)
802 {
803 if (this->headersAreRead)
804 {
805 int numberOfLevels = this->GetNumberOfLevels() + 1;
806 int levelBlocksLo = 0;
807 int levelBlocksHi = 0;
808 for (int cc = 0; cc < numberOfLevels; ++cc)
809 {
810 levelBlocksHi += this->LevelHeader[cc]->levelBoxArraySize;
811 if (blockIdx >= levelBlocksLo && blockIdx < levelBlocksHi)
812 {
813 return cc;
814 }
815 levelBlocksLo = levelBlocksHi;
816 }
817 }
818 return (-1);
819 }
820
GetNumberOfBlocks()821 int vtkAMReXGridReaderInternal::GetNumberOfBlocks()
822 {
823 if (this->headersAreRead)
824 {
825 int numberOfLevels = this->GetNumberOfLevels() + 1;
826 int numberOfBlocks = 0;
827 for (int i = 0; i < numberOfLevels; ++i)
828 {
829 numberOfBlocks += this->Header->levelSize[i];
830 }
831 return numberOfBlocks;
832 }
833 return (-1);
834 }
835
GetBlockIndexWithinLevel(int blockIdx,int level)836 int vtkAMReXGridReaderInternal::GetBlockIndexWithinLevel(int blockIdx, int level)
837 {
838 if (this->headersAreRead)
839 {
840 int blockIndexWithinLevel = blockIdx;
841 for (int i = 0; i < level; ++i)
842 {
843 blockIndexWithinLevel -= this->Header->levelSize[i];
844 }
845 return blockIndexWithinLevel;
846 }
847 return (-1);
848 }
849
GetBlockAttribute(const char * attribute,int blockIdx,vtkDataSet * pDataSet)850 void vtkAMReXGridReaderInternal::GetBlockAttribute(
851 const char* attribute, int blockIdx, vtkDataSet* pDataSet)
852 {
853 if (this->headersAreRead)
854 {
855 if (attribute == nullptr || blockIdx < 0 || pDataSet == nullptr ||
856 blockIdx >= this->GetNumberOfBlocks())
857 {
858 return;
859 }
860
861 //
862 // orders.
863 //
864 // int big_float_order[] = { 1, 2, 3, 4 };
865 int little_float_order[] = { 4, 3, 2, 1 };
866 // int mid_float_order_2[] = { 2, 1, 4, 3 };
867 // int big_double_order[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
868 int little_double_order[] = { 8, 7, 6, 5, 4, 3, 2, 1 };
869 // int mid_double_order_2[] = { 2, 1, 4, 3, 6, 5, 8, 7 };
870 //
871 // formats.
872 //
873 long ieee_float[] = { 32L, 8L, 23L, 0L, 1L, 9L, 0L, 0x7FL };
874 long ieee_double[] = { 64L, 11L, 52L, 0L, 1L, 12L, 0L, 0x3FFL };
875
876 int offsetOfAttribute = this->GetOffsetOfAttribute(attribute);
877 int theLevel = this->GetBlockLevel(blockIdx);
878 int blockIdxWithinLevel = this->GetBlockIndexWithinLevel(blockIdx, theLevel);
879 if (debugReader)
880 std::cout << "blockIdx " << blockIdx << " attribute " << attribute << " offset of attribute "
881 << offsetOfAttribute << " Level " << theLevel << " blockIdx within Level "
882 << blockIdxWithinLevel << std::endl;
883 const std::string FABFileName = this->FileName + "/" + this->Header->levelPrefix[theLevel] +
884 "/" + this->LevelHeader[theLevel]->levelFABFile[blockIdxWithinLevel];
885 if (debugReader)
886 std::cout << "FABFile " << FABFileName << " Offset "
887 << this->LevelHeader[theLevel]->levelFileOffset[blockIdxWithinLevel] << std::endl;
888 std::filebuf fb;
889 if (fb.open(FABFileName, std::ios::binary | std::ios::in))
890 {
891 std::istream is(&fb);
892 is.seekg(this->LevelHeader[theLevel]->levelFileOffset[blockIdxWithinLevel]);
893 //
894 // Read FAB Header
895 //
896 this->ReadFAB(is);
897 // int version =
898 this->ReadVersion(is);
899 int dimension = this->Header->dim;
900 RealDescriptor* ird = this->ReadRealDescriptor(is);
901 std::vector<int> boxArray(3 * dimension);
902 std::vector<int> boxArrayDim(dimension);
903 int numberOfPoints = ReadBoxArray(is, boxArray.data(), boxArrayDim.data());
904 // int numberOfAttributes =
905 this->ReadNumberOfAttributes(is);
906
907 //
908 // Skip the Line Feed (linefeed+1)
909 // Jump to the desired attribute (offsetOfAttribute*(numberOfPoints*ird->numBytes()))
910 // - Patrick O'Leary
911 //
912 int linefeed = is.tellg();
913
914 if (debugReader)
915 {
916 for (int i = 0; i < dimension; ++i)
917 std::cout << boxArrayDim[i] << " ";
918 std::cout << std::endl;
919 }
920
921 // read every component of the variable into the buffers vector
922 std::string attributeName(attribute);
923 int nComps = static_cast<int>(this->Header->parsedVariableNames[attributeName].size());
924 std::vector<std::vector<char>> buffers(nComps);
925 for (int i = 0; i < nComps; ++i)
926 {
927 int compIndex = this->Header->parsedVariableNames[attributeName][i];
928 std::string compName = this->Header->variableNames[compIndex];
929 offsetOfAttribute = this->GetOffsetOfAttribute(compName.c_str());
930 is.seekg((linefeed + 1) + (offsetOfAttribute * (numberOfPoints * ird->numBytes())));
931 buffers[i].resize(numberOfPoints * ird->numBytes());
932 this->ReadBlockAttribute(is, numberOfPoints, ird->numBytes(), buffers[i].data());
933 }
934
935 RealDescriptor* ord = nullptr;
936 // copy buffers into vtkAOSDataArrayTemplate
937 if (ird->numBytes() == 4)
938 {
939 vtkNew<vtkAOSDataArrayTemplate<float>> dataArray;
940 ord = new RealDescriptor(ieee_float, little_float_order, 4);
941 this->CreateVTKAttributeArray(
942 dataArray.Get(), ord, ird, buffers, numberOfPoints, attributeName);
943 pDataSet->GetCellData()->AddArray(dataArray);
944 }
945 else
946 {
947 vtkNew<vtkAOSDataArrayTemplate<double>> dataArray;
948 ord = new RealDescriptor(ieee_double, little_double_order, 8);
949 this->CreateVTKAttributeArray(
950 dataArray.Get(), ord, ird, buffers, numberOfPoints, attributeName);
951 pDataSet->GetCellData()->AddArray(dataArray);
952 }
953 delete ord;
954 delete ird;
955
956 if (debugReader)
957 {
958 std::cout << is.tellg() << " "
959 << this->LevelHeader[theLevel]->levelFileOffset[blockIdxWithinLevel] << " "
960 << numberOfPoints << std::endl;
961 }
962 fb.close();
963 }
964 }
965 }
966
GetOffsetOfAttribute(const char * attribute)967 int vtkAMReXGridReaderInternal::GetOffsetOfAttribute(const char* attribute)
968 {
969 int i = 0, position = 0;
970 bool found = false;
971
972 while (i < this->Header->variableNamesSize && !found)
973 {
974 if (strcmp(this->Header->variableNames[i].c_str(), attribute) == 0)
975 {
976 found = true;
977 position = i;
978 }
979 ++i;
980 }
981 if (found)
982 {
983 return position;
984 }
985 else
986 return (-1);
987 }
988
ReadFAB(std::istream & is)989 void vtkAMReXGridReaderInternal::ReadFAB(std::istream& is)
990 {
991 char f, a, b;
992 is >> f;
993 is >> a;
994 is >> b;
995 if (debugReader)
996 std::cout << f << a << b;
997 }
998
ReadVersion(std::istream & is)999 int vtkAMReXGridReaderInternal::ReadVersion(std::istream& is)
1000 {
1001 char colon;
1002 is >> colon;
1003 if (colon == ':')
1004 {
1005 if (debugReader)
1006 std::cout << colon << "!" << std::endl;
1007
1008 return 0;
1009 }
1010 else
1011 {
1012 is.putback(colon);
1013 if (debugReader)
1014 std::cout << " ";
1015 return 1;
1016 }
1017 }
1018
ReadOrder(std::istream & is,std::vector<int> & ar)1019 void vtkAMReXGridReaderInternal::ReadOrder(std::istream& is, std::vector<int>& ar)
1020 {
1021 char c;
1022 is >> c; // '('
1023 int size;
1024 is >> size;
1025 is >> c; // ','
1026 is >> c; // '('
1027 ar.resize(size);
1028 for (int i = 0; i < size; ++i)
1029 is >> ar[i];
1030 is >> c; // ')'
1031 is >> c; // ')'
1032 }
1033
PrintOrder(std::vector<int> & ar)1034 void vtkAMReXGridReaderInternal::PrintOrder(std::vector<int>& ar)
1035 {
1036 size_t size = ar.size();
1037 std::cout << "(" << size << ", (";
1038 for (size_t i = 0; i < size; ++i)
1039 {
1040 std::cout << ar[i];
1041 if (i < size - 1)
1042 std::cout << " ";
1043 }
1044 std::cout << "))";
1045 }
1046
ReadFormat(std::istream & is,std::vector<long> & ar)1047 void vtkAMReXGridReaderInternal::ReadFormat(std::istream& is, std::vector<long>& ar)
1048 {
1049 char c;
1050 is >> c; // '('
1051 int size;
1052 is >> size;
1053 is >> c; // ','
1054 is >> c; // '('
1055 ar.resize(size);
1056 for (int i = 0; i < size; ++i)
1057 is >> ar[i];
1058 is >> c; // ')'
1059 is >> c; // ')'
1060 }
1061
PrintFormat(std::vector<long> & ar)1062 void vtkAMReXGridReaderInternal::PrintFormat(std::vector<long>& ar)
1063 {
1064 size_t size = ar.size();
1065 std::cout << "(" << size << ", (";
1066 for (size_t i = 0; i < size; ++i)
1067 {
1068 std::cout << ar[i];
1069 if (i < size - 1)
1070 std::cout << " ";
1071 }
1072 std::cout << "))";
1073 }
1074
ReadRealDescriptor(std::istream & is)1075 RealDescriptor* vtkAMReXGridReaderInternal::ReadRealDescriptor(std::istream& is)
1076 {
1077 std::vector<long> fmt;
1078 std::vector<int> ord;
1079 char c;
1080 is >> c; // '('
1081 if (debugReader)
1082 std::cout << c;
1083 this->ReadFormat(is, fmt);
1084 if (debugReader)
1085 this->PrintFormat(fmt);
1086 is >> c; // ','
1087 if (debugReader)
1088 std::cout << c;
1089 this->ReadOrder(is, ord);
1090 if (debugReader)
1091 this->PrintOrder(ord);
1092 is >> c; // ')'
1093 if (debugReader)
1094 std::cout << c;
1095 //
1096 // ord.size() is either 4 or 8 for float or double respectively - cast to int is safe
1097 //
1098 return new RealDescriptor(&fmt[0], &ord[0], static_cast<int>(ord.size()));
1099 }
1100
ReadBoxArray(std::istream & is,int * boxArray,int * boxArrayDim)1101 int vtkAMReXGridReaderInternal::ReadBoxArray(std::istream& is, int* boxArray, int* boxArrayDim)
1102 {
1103 char c;
1104 is >> c; // read '('
1105 for (int dd = 0; dd < 3; ++dd)
1106 {
1107 is >> c; // read '('
1108 for (int ee = 0; ee < this->Header->dim; ++ee)
1109 {
1110 is >> boxArray[this->Header->dim * dd + ee];
1111 if (ee < (this->Header->dim - 1))
1112 {
1113 is >> c; // read ','
1114 }
1115 else
1116 {
1117 is >> c; // read ')'
1118 }
1119 }
1120 }
1121 is >> c; // read ')'
1122
1123 //
1124 // block dimension - '(hi - lo + 1)' is the number of cells '+ 1' is the number of points
1125 //
1126 int numberOfPoints = 1;
1127 for (int i = 0; i < this->Header->dim; ++i)
1128 {
1129 boxArrayDim[i] =
1130 ((boxArray[this->Header->dim * 1 + i] - boxArray[this->Header->dim * 0 + i]) + 1);
1131 numberOfPoints *= boxArrayDim[i];
1132 }
1133 if (debugReader)
1134 this->PrintBoxArray(boxArray);
1135 return numberOfPoints;
1136 }
1137
PrintBoxArray(int * boxArray)1138 void vtkAMReXGridReaderInternal::PrintBoxArray(int* boxArray)
1139 {
1140 std::cout << "(";
1141 std::cout << "(";
1142 for (int space = 0; space < this->Header->dim; ++space)
1143 {
1144 std::cout << boxArray[0 * this->Header->dim + space];
1145 if (space < (this->Header->dim - 1))
1146 std::cout << ",";
1147 }
1148 std::cout << ") ";
1149 std::cout << "(";
1150 for (int space = 0; space < this->Header->dim; ++space)
1151 {
1152 std::cout << boxArray[1 * this->Header->dim + space];
1153 if (space < (this->Header->dim - 1))
1154 std::cout << ",";
1155 }
1156 std::cout << ") ";
1157 std::cout << "(";
1158 for (int space = 0; space < this->Header->dim; ++space)
1159 {
1160 std::cout << boxArray[2 * this->Header->dim + space];
1161 if (space < (this->Header->dim - 1))
1162 std::cout << ",";
1163 }
1164 std::cout << ")";
1165 std::cout << ")";
1166 }
1167
ReadNumberOfAttributes(std::istream & is)1168 int vtkAMReXGridReaderInternal::ReadNumberOfAttributes(std::istream& is)
1169 {
1170 int numberOfAttributes;
1171 is >> numberOfAttributes;
1172 if (debugReader)
1173 std::cout << " " << numberOfAttributes << std::endl;
1174 return numberOfAttributes;
1175 }
1176
ReadBlockAttribute(std::istream & is,int numberOfPoints,int size,char * buffer)1177 void vtkAMReXGridReaderInternal::ReadBlockAttribute(
1178 std::istream& is, int numberOfPoints, int size, char* buffer)
1179 {
1180 is.read(buffer, numberOfPoints * size);
1181 }
1182
Convert(void * out,const void * in,long nitems,const RealDescriptor & ord,const RealDescriptor & ird)1183 void vtkAMReXGridReaderInternal::Convert(
1184 void* out, const void* in, long nitems, const RealDescriptor& ord, const RealDescriptor& ird)
1185 {
1186 if (ord == ird)
1187 {
1188 size_t n = size_t(nitems);
1189 memcpy(out, in, n * ord.numBytes());
1190 }
1191 else if (ord.formatarray() == ird.formatarray())
1192 {
1193 this->PermuteOrder(out, in, nitems, ord.order(), ird.order(), ord.numBytes());
1194 }
1195 else
1196 {
1197 //
1198 // We don't handle changing real formats
1199 //
1200 }
1201 }
1202
PermuteOrder(void * out,const void * in,long nitems,const int * outord,const int * inord,int REALSIZE)1203 void vtkAMReXGridReaderInternal::PermuteOrder(
1204 void* out, const void* in, long nitems, const int* outord, const int* inord, int REALSIZE)
1205 {
1206 char const* pin = static_cast<char const*>(in);
1207 char* pout = static_cast<char*>(out);
1208
1209 pin--;
1210 pout--;
1211
1212 for (; nitems > 0; nitems--, pin += REALSIZE, pout += REALSIZE)
1213 {
1214 for (int i = 0; i < REALSIZE; i++)
1215 pout[outord[i]] = pin[inord[i]];
1216 }
1217 }
1218