1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMFIXReader.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 // Thanks to Phil Nicoletti, Terry Jordan and Brian Dotson at the
16 // National Energy Technology Laboratory who developed this class.
17 // Please address all comments to Terry Jordan (terry.jordan@netl.doe.gov)
18 //
19
20 #include "vtkMFIXReader.h"
21
22 #include "vtkInformation.h"
23 #include "vtkInformationVector.h"
24 #include "vtkObjectFactory.h"
25 #include "vtkErrorCode.h"
26 #include "vtkUnstructuredGrid.h"
27 #include "vtkPointData.h"
28 #include "vtkCellData.h"
29 #include "vtkDoubleArray.h"
30 #include "vtkIntArray.h"
31 #include "vtkCellArray.h"
32 #include "vtkHexahedron.h"
33 #include "vtkFloatArray.h"
34 #include <string>
35 #include "vtkDataArraySelection.h"
36 #include "vtkWedge.h"
37 #include "vtkQuad.h"
38 #include "vtkStreamingDemandDrivenPipeline.h"
39 #include "vtkStringArray.h"
40
41 vtkStandardNewMacro(vtkMFIXReader);
42
43 //----------------------------------------------------------------------------
vtkMFIXReader()44 vtkMFIXReader::vtkMFIXReader()
45 {
46 this->FileName = NULL;
47 this->NumberOfCells = 0;
48 this->NumberOfPoints = 0;
49 this->NumberOfCellFields = 0;
50 this->RequestInformationFlag = 0;
51 this->MakeMeshFlag = 0;
52 this->Minimum = vtkFloatArray::New();
53 this->Maximum = vtkFloatArray::New();
54 this->VectorLength = vtkIntArray::New();
55 this->CellDataArray = NULL;
56 this->DimensionIc = 5;
57 this->DimensionBc = 5;
58 this->DimensionC = 5;
59 this->DimensionIs = 5;
60 this->NumberOfSPXFilesUsed = 9;
61 this->NumberOfScalars = 0;
62 this->BkEpsilon = false;
63 this->NumberOfReactionRates = 0;
64 this->FileExtension[0] = '1';
65 this->FileExtension[1] = '2';
66 this->FileExtension[2] = '3';
67 this->FileExtension[3] = '4';
68 this->FileExtension[4] = '5';
69 this->FileExtension[5] = '6';
70 this->FileExtension[6] = '7';
71 this->FileExtension[7] = '8';
72 this->FileExtension[8] = '9';
73 this->FileExtension[9] = 'A';
74 this->FileExtension[10] = 'B';
75 this->VersionNumber = 0;
76
77 this->CellDataArraySelection = vtkDataArraySelection::New();
78 this->Points = vtkPoints::New();
79 this->Mesh = vtkUnstructuredGrid::New();
80 this->AHexahedron = vtkHexahedron::New();
81 this->AQuad = vtkQuad::New();
82 this->AWedge = vtkWedge::New();
83 this->NMax = vtkIntArray::New();
84 this->C = vtkDoubleArray::New();
85 this->Dx = vtkDoubleArray::New();
86 this->Dy = vtkDoubleArray::New();
87 this->Dz = vtkDoubleArray::New();
88 this->TempI = vtkIntArray::New();
89 this->TempD = vtkDoubleArray::New();
90 this->Flag = vtkIntArray::New();
91 this->VariableNames = vtkStringArray::New();
92 this->VariableComponents = vtkIntArray::New();
93 this->VariableIndexToSPX = vtkIntArray::New();
94 this->VariableTimesteps = vtkIntArray::New();
95 this->VariableTimestepTable = vtkIntArray::New();
96 this->SPXToNVarTable = vtkIntArray::New();
97 this->VariableToSkipTable = vtkIntArray::New();
98 this->SpxFileExists = vtkIntArray::New();
99 this->SetNumberOfInputPorts(0);
100 this->SPXTimestepIndexTable = vtkIntArray::New();
101
102 // Time support:
103 this->TimeStep = 0; // By default the file does not have timestep
104 this->TimeStepRange[0] = 0;
105 this->TimeStepRange[1] = 0;
106 this->NumberOfTimeSteps = 1;
107 this->TimeSteps = 0;
108 this->CurrentTimeStep = 0;
109 this->TimeStepWasReadOnce = 0;
110 }
111
112 //----------------------------------------------------------------------------
~vtkMFIXReader()113 vtkMFIXReader::~vtkMFIXReader()
114 {
115 delete [] this->FileName;
116
117 if( this->CellDataArray )
118 {
119 for (int j = 0; j <= this->VariableNames->GetMaxId(); j++)
120 {
121 this->CellDataArray[j]->Delete();
122 }
123 delete [] this->CellDataArray;
124 }
125
126 this->CellDataArraySelection->Delete();
127 this->Points->Delete();
128 this->Mesh->Delete();
129 this->AHexahedron->Delete();
130 this->AWedge->Delete();
131 this->AQuad->Delete();
132 this->NMax->Delete();
133 this->C->Delete();
134 this->Dx->Delete();
135 this->Dy->Delete();
136 this->Dz->Delete();
137 this->TempI->Delete();
138 this->TempD->Delete();
139 this->Flag->Delete();
140 this->VariableNames->Delete();
141 this->VariableComponents->Delete();
142 this->VariableIndexToSPX->Delete();
143 this->VariableTimesteps->Delete();
144 this->VariableTimestepTable->Delete();
145 this->SPXToNVarTable->Delete();
146 this->VariableToSkipTable->Delete();
147 this->SpxFileExists->Delete();
148 this->Minimum->Delete();
149 this->Maximum->Delete();
150 this->VectorLength->Delete();
151 this->SPXTimestepIndexTable->Delete();
152 }
153
154 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)155 int vtkMFIXReader::RequestData(
156 vtkInformation *vtkNotUsed(request),
157 vtkInformationVector **vtkNotUsed(inputVector),
158 vtkInformationVector *outputVector)
159 {
160 vtkInformation* outInfo = outputVector->GetInformationObject(0);
161 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
162 outInfo->Get(vtkDataObject::DATA_OBJECT()));
163 vtkDebugMacro( << "Reading MFIX file");
164
165 // Save the time value in the output data information.
166 int length = outInfo->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
167 double* steps = outInfo->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
168
169 if (outInfo->Has( vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()))
170 {
171 // Get the requested time step. We only support requests of a single time
172 // step in this reader right now
173 double requestedTimeStep =
174 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP());
175
176
177 //find the timestep with the closest value
178 int cnt=0;
179 int closestStep=0;
180 double minDist=-1;
181 for (cnt=0;cnt<length;cnt++)
182 {
183 double tdist=(steps[cnt]-requestedTimeStep>requestedTimeStep-steps[cnt])?
184 steps[cnt]-requestedTimeStep:
185 requestedTimeStep-steps[cnt];
186 if (minDist<0 || tdist<minDist)
187 {
188 minDist=tdist;
189 closestStep=cnt;
190 }
191 }
192 this->CurrentTimeStep=closestStep;
193 }
194 else
195 {
196 this->CurrentTimeStep = this->TimeStep;
197 }
198
199 this->MakeMesh(output);
200 output->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(), steps[this->CurrentTimeStep]);
201 return 1;
202 }
203
204 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)205 void vtkMFIXReader::PrintSelf(ostream& os, vtkIndent indent)
206 {
207 this->Superclass::PrintSelf(os,indent);
208 os << indent << "File Name: "
209 << (this->FileName ? this->FileName : "(none)") << "\n";
210 os << indent << "Number Of Nodes: " << this->NumberOfPoints << endl;
211 os << indent << "Number Of Cells: " << this->NumberOfCells << endl;
212 os << indent << "Number Of Cell Fields: " << this->NumberOfCellFields << endl;
213 os << indent << "Time Step Range: "
214 << this->TimeStepRange[0] << " - " << this->TimeStepRange[1]
215 << endl;
216 os << indent << "Time Step: " << this->TimeStep << endl;
217 os << indent << "Number of Time Steps: " << this->NumberOfTimeSteps << endl;
218 }
219
220 //----------------------------------------------------------------------------
MakeMesh(vtkUnstructuredGrid * output)221 void vtkMFIXReader::MakeMesh(vtkUnstructuredGrid *output)
222 {
223 output->Allocate();
224
225 if (this->MakeMeshFlag == 0)
226 {
227 //Points->SetNumberOfPoints((this->IMaximum2+1)
228 // *(this->JMaximum2+1)*(this->KMaximum2+1));
229
230 //
231 // Cartesian type mesh
232 //
233 if ( !strcmp(this->CoordinateSystem,"CARTESIAN") && (this->KMaximum2 != 1))
234 {
235 double pointX = -this->Dx->GetValue(0);
236 double pointY = -this->Dy->GetValue(0);
237 double pointZ = -this->Dz->GetValue(0);
238 for (int k = 0; k <= this->KMaximum2; k++)
239 {
240 for (int j = 0; j <= this->JMaximum2; j++)
241 {
242 for (int i = 0; i <= this->IMaximum2; i++)
243 {
244 this->Points->InsertNextPoint( pointX, pointY, pointZ );
245 if ( i == this->IMaximum2 )
246 {
247 pointX = pointX + this->Dx->GetValue(i-1);
248 }
249 else
250 {
251 pointX = pointX + this->Dx->GetValue(i);
252 }
253 }
254 pointX = -this->Dx->GetValue(0);
255 if ( j == this->JMaximum2)
256 {
257 pointY = pointY + this->Dy->GetValue(j-1);
258 }
259 else
260 {
261 pointY = pointY + this->Dy->GetValue(j);
262 }
263 }
264 pointY = -this->Dy->GetValue(0);
265 if ( k == this->KMaximum2)
266 {
267 pointZ = pointZ + this->Dz->GetValue(k-1);
268 }
269 else
270 {
271 pointZ = pointZ + this->Dz->GetValue(k);
272 }
273 }
274 }
275 else if (!strcmp(this->CoordinateSystem,"CARTESIAN") &&
276 (this->KMaximum2 == 1))
277 {
278 double pointX = -this->Dx->GetValue(0);
279 double pointY = -this->Dy->GetValue(0);
280 double pointZ = 0.0;
281 for (int j = 0; j <= this->JMaximum2; j++)
282 {
283 for (int i = 0; i <= this->IMaximum2; i++)
284 {
285 this->Points->InsertNextPoint(pointX, pointY, pointZ );
286 if ( i == this->IMaximum2 )
287 {
288 pointX = pointX + this->Dx->GetValue(i-1);
289 }
290 else
291 {
292 pointX = pointX + this->Dx->GetValue(i);
293 }
294 }
295 pointX = -this->Dx->GetValue(0);
296 if ( j == this->JMaximum2)
297 {
298 pointY = pointY + this->Dy->GetValue(j-1);
299 }
300 else
301 {
302 pointY = pointY + this->Dy->GetValue(j);
303 }
304 }
305 }
306 else if (!strcmp(this->CoordinateSystem,"CYLINDRICAL") &&
307 (this->KMaximum2 == 1))
308 {
309 double pointX = -this->Dx->GetValue(0);
310 double pointY = -this->Dy->GetValue(0);
311 double pointZ = 0.0;
312 for (int j = 0; j <= this->JMaximum2; j++)
313 {
314 for (int i = 0; i <= this->IMaximum2; i++)
315 {
316 this->Points->InsertNextPoint( pointX, pointY, pointZ );
317 if ( i == this->IMaximum2 )
318 {
319 pointX = pointX + this->Dx->GetValue(i-1);
320 }
321 else if ( i == 0 )
322 {
323 pointX = 0;
324 }
325 else
326 {
327 pointX = pointX + this->Dx->GetValue(i);
328 }
329 }
330 pointX = -this->Dx->GetValue(0);
331 if ( j == this->JMaximum2)
332 {
333 pointY = pointY + this->Dy->GetValue(j-1);
334 }
335 else
336 {
337 pointY = pointY + this->Dy->GetValue(j);
338 }
339 }
340 }
341 else
342 {
343 //
344 // Cylindrical Type Mesh
345 //
346 double pointX = -this->Dx->GetValue(0);
347 double pointY = -this->Dy->GetValue(0);
348 double pointZ = -this->Dz->GetValue(0);
349 double radialX = 0.0;
350 double radialY = 0.0;
351 double radialZ = 0.0;
352 for (int k = 0; k <= this->KMaximum2; k++)
353 {
354 for (int j = 0; j <= this->JMaximum2; j++)
355 {
356 for (int i = 0; i <= this->IMaximum2; i++)
357 {
358 this->Points->InsertNextPoint( radialX, radialY, radialZ );
359 if ( i == this->IMaximum2 )
360 {
361 pointX = pointX + this->Dx->GetValue(i-1);
362 }
363 else if ( i == 0 )
364 {
365 pointX = -this->Dx->GetValue(0);
366 }
367 else
368 {
369 pointX = pointX + this->Dx->GetValue(i);
370 }
371 radialX = pointX * cos(pointZ);
372 radialZ = pointX * sin(pointZ) * -1;
373 }
374 pointX = -this->Dx->GetValue(0);
375 radialX = 0.0;
376 radialZ = 0.0;
377 if ( j == this->JMaximum2)
378 {
379 pointY = pointY + this->Dy->GetValue(j-1);
380 }
381 else
382 {
383 pointY = pointY + this->Dy->GetValue(j);
384 }
385 radialY = pointY;
386 }
387 pointY = -this->Dy->GetValue(0);
388 radialY = 0.0;
389 if ( k == this->KMaximum2)
390 {
391 pointZ = pointZ + this->Dz->GetValue(k-1);
392 }
393 else
394 {
395 pointZ = pointZ + this->Dz->GetValue(k);
396 }
397 }
398 }
399
400 //
401 // Let's put the points in a mesh
402 //
403 this->Mesh->SetPoints(this->Points);
404 int point0 = 0;
405 int count = 0;
406 if (!strcmp(this->CoordinateSystem,"CYLINDRICAL") &&
407 (this->KMaximum2 == 1))
408 {
409 for (int j = 0; j < this->JMaximum2; j++)
410 {
411 for (int i = 0; i < this->IMaximum2; i++)
412 {
413 if ( this->Flag->GetValue(count) < 10 )
414 {
415 this->AQuad->GetPointIds()->SetId( 0, point0);
416 this->AQuad->GetPointIds()->SetId( 1, point0+1);
417 this->AQuad->GetPointIds()->SetId( 2, point0+2+this->IMaximum2);
418 this->AQuad->GetPointIds()->SetId( 3, point0+1+this->IMaximum2);
419 this->Mesh->InsertNextCell(this->AQuad->GetCellType(),
420 this->AQuad->GetPointIds());
421 }
422 point0++;
423 count++;
424 }
425 point0++;
426 }
427 point0 = point0 + this->IMaximum2+1;
428 }
429 else if (!strcmp(this->CoordinateSystem,"CYLINDRICAL") &&
430 (this->KMaximum2 != 1))
431 {
432 for (int k = 0; k < this->KMaximum2; k++)
433 {
434 for (int j = 0; j < this->JMaximum2; j++)
435 {
436 for (int i = 0; i < this->IMaximum2; i++)
437 {
438 if ( this->Flag->GetValue(count) < 10 )
439 {
440 if (( k == (this->KMaximum2-2)) && (i != 1))
441 {
442 this->AHexahedron->GetPointIds()->SetId( 0, point0);
443 this->AHexahedron->GetPointIds()->SetId( 1, point0+1);
444 this->AHexahedron->GetPointIds()->SetId( 2,
445 (point0+1+((this->IMaximum2+1)*(this->JMaximum2+1)))-
446 ((this->IMaximum2+1)*(this->JMaximum2+1)
447 *(this->KMaximum2-2)));
448 this->AHexahedron->GetPointIds()->SetId( 3,
449 (point0+((this->IMaximum2+1)*(this->JMaximum2+1)))-
450 ((this->IMaximum2+1)*(this->JMaximum2+1)
451 *(this->KMaximum2-2)));
452 this->AHexahedron->GetPointIds()->
453 SetId( 4, point0+1+this->IMaximum2);
454 this->AHexahedron->GetPointIds()->
455 SetId( 5, point0+2+this->IMaximum2);
456 this->AHexahedron->GetPointIds()->
457 SetId( 6, (point0+2+this->IMaximum2 +
458 ((this->IMaximum2+1)*(this->JMaximum2+1))) -
459 ((this->IMaximum2+1)*(this->JMaximum2+1)
460 * (this->KMaximum2-2)));
461 this->AHexahedron->GetPointIds()->
462 SetId( 7, (point0+1+this->IMaximum2 +
463 ((this->IMaximum2+1)*(this->JMaximum2+1)))-
464 ((this->IMaximum2+1)*(this->JMaximum2+1)
465 *(this->KMaximum2-2)));
466 this->Mesh->InsertNextCell(this->AHexahedron->GetCellType(),
467 this->AHexahedron->GetPointIds());
468 }
469 else if ((k != (this->KMaximum2-2)) && (i != 1))
470 {
471 this->AHexahedron->GetPointIds()->SetId( 0, point0);
472 this->AHexahedron->GetPointIds()->SetId( 1, point0+1);
473 this->AHexahedron->GetPointIds()->SetId( 2,
474 point0+1+((this->IMaximum2+1)*(this->JMaximum2+1)));
475 this->AHexahedron->GetPointIds()->SetId( 3,
476 point0+((this->IMaximum2+1)*(this->JMaximum2+1)));
477 this->AHexahedron->GetPointIds()->
478 SetId( 4, point0+1+this->IMaximum2);
479 this->AHexahedron->GetPointIds()->
480 SetId( 5, point0+2+this->IMaximum2);
481 this->AHexahedron->GetPointIds()->
482 SetId( 6, point0+2+this->IMaximum2+
483 ((this->IMaximum2+1)*(this->JMaximum2+1)));
484 this->AHexahedron->GetPointIds()->
485 SetId( 7, point0+1+this->IMaximum2+
486 ((this->IMaximum2+1)*(this->JMaximum2+1)));
487 this->Mesh->InsertNextCell(this->AHexahedron->GetCellType(),
488 this->AHexahedron->GetPointIds());
489 }
490 else if ( (k != (this->KMaximum2-2)) && (i == 1))
491 {
492 this->AWedge->GetPointIds()->SetId( 0, j*(this->IMaximum2+1));
493 this->AWedge->GetPointIds()->SetId( 1, point0+1);
494 this->AWedge->GetPointIds()->
495 SetId( 2, point0+1+((this->IMaximum2+1)
496 *(this->JMaximum2+1)));
497 this->AWedge->GetPointIds()->
498 SetId( 3, (j+1)*(this->IMaximum2+1));
499 this->AWedge->GetPointIds()->
500 SetId( 4, point0+2+this->IMaximum2);
501 this->AWedge->GetPointIds()->
502 SetId( 5, point0+2+this->IMaximum2+
503 ((this->IMaximum2+1)*(this->JMaximum2+1)));
504 this->Mesh->InsertNextCell(this->AWedge->GetCellType(),
505 this->AWedge->GetPointIds());
506 }
507 else if (( k == (this->KMaximum2-2)) && (i == 1))
508 {
509 this->AWedge->GetPointIds()->SetId( 0, j*(this->IMaximum2+1));
510 this->AWedge->GetPointIds()->SetId( 1, point0+1);
511 this->AWedge->GetPointIds()->SetId( 2,
512 (point0+1+((this->IMaximum2+1)
513 *(this->JMaximum2+1)))-((this->IMaximum2+1)
514 *(this->JMaximum2+1)*(this->KMaximum2-2)));
515 this->AWedge->GetPointIds()->
516 SetId( 3, (j+1)*(this->IMaximum2+1));
517 this->AWedge->GetPointIds()->
518 SetId( 4, point0+2+this->IMaximum2);
519 this->AWedge->GetPointIds()->
520 SetId( 5, (point0+2+this->IMaximum2 +
521 ((this->IMaximum2+1)*(this->JMaximum2+1)))
522 -((this->IMaximum2+1)
523 *(this->JMaximum2+1)*(this->KMaximum2-2)));
524 this->Mesh->InsertNextCell(this->AWedge->GetCellType(),
525 this->AWedge->GetPointIds());
526 }
527 }
528 point0++;
529 count++;
530 }
531 point0++;
532 }
533 point0 = point0 + this->IMaximum2+1;
534 }
535 }
536 else if (!strcmp(this->CoordinateSystem,"CARTESIAN") &&
537 (this->KMaximum2 == 1))
538 {
539 for (int j = 0; j < this->JMaximum2; j++)
540 {
541 for (int i = 0; i < this->IMaximum2; i++)
542 {
543 if ( this->Flag->GetValue(count) < 10 )
544 {
545 this->AQuad->GetPointIds()->SetId( 0, point0);
546 this->AQuad->GetPointIds()->SetId( 1, point0+1);
547 this->AQuad->GetPointIds()->SetId( 2, point0+2+this->IMaximum2);
548 this->AQuad->GetPointIds()->SetId( 3, point0+1+this->IMaximum2);
549 this->Mesh->InsertNextCell(this->AQuad->GetCellType(),
550 this->AQuad->GetPointIds());
551 }
552 point0++;
553 count++;
554 }
555 point0++;
556 }
557 }
558 else
559 {
560 for (int k = 0; k < this->KMaximum2; k++)
561 {
562 for (int j = 0; j < this->JMaximum2; j++)
563 {
564 for (int i = 0; i < this->IMaximum2; i++)
565 {
566 if ( this->Flag->GetValue(count) < 10 )
567 {
568 this->AHexahedron->GetPointIds()->SetId( 0, point0);
569 this->AHexahedron->GetPointIds()->SetId( 1, point0+1);
570 this->AHexahedron->GetPointIds()->
571 SetId( 2, point0+1+((this->IMaximum2+1)
572 *(this->JMaximum2+1)));
573 this->AHexahedron->GetPointIds()->
574 SetId( 3, point0+((this->IMaximum2+1)
575 *(this->JMaximum2+1)));
576 this->AHexahedron->GetPointIds()->
577 SetId( 4, point0+1+this->IMaximum2);
578 this->AHexahedron->GetPointIds()->
579 SetId( 5, point0+2+this->IMaximum2);
580 this->AHexahedron->GetPointIds()->
581 SetId( 6, point0+2+this->IMaximum2 +
582 ((this->IMaximum2+1)*(this->JMaximum2+1)));
583 this->AHexahedron->GetPointIds()->
584 SetId( 7, point0+1+this->IMaximum2 +
585 ((this->IMaximum2+1)*(this->JMaximum2+1)));
586 this->Mesh->InsertNextCell(this->AHexahedron->GetCellType(),
587 this->AHexahedron->GetPointIds());
588 }
589 point0++;
590 count++;
591 }
592 point0++;
593 }
594 point0 = point0 + this->IMaximum2+1;
595 }
596 }
597
598
599 this->CellDataArray = new vtkFloatArray
600 * [this->VariableNames->GetMaxId()+2];
601 for (int j = 0; j <= this->VariableNames->GetMaxId(); j++)
602 {
603 this->CellDataArray[ j ] = vtkFloatArray::New();
604 this->CellDataArray[ j ]->SetName(this->VariableNames->GetValue(j));
605 this->CellDataArray[ j ]->
606 SetNumberOfComponents(this->VariableComponents->GetValue(j));
607 }
608
609 this->MakeMeshFlag = 1;
610 }
611
612 output->DeepCopy(this->Mesh); // If mesh has already
613 // been made copy it to output
614 int first = 0;
615 for (int j = 0; j <= this->VariableNames->GetMaxId(); j++)
616 {
617 if ( this->CellDataArraySelection->GetArraySetting(j) == 1 )
618 {
619 if (this->VariableComponents->GetValue(j) == 1)
620 {
621 this->GetVariableAtTimestep( j, this->CurrentTimeStep, CellDataArray[j]);
622 }
623 else
624 {
625 if ( !strcmp(CoordinateSystem,"CYLINDRICAL" ))
626 {
627 this->ConvertVectorFromCylindricalToCartesian( j-3, j-1);
628 }
629 this->FillVectorVariable( j-3, j-2, j-1, this->CellDataArray[j]);
630 }
631 if (first == 0)
632 {
633 output->GetCellData()->SetScalars(this->CellDataArray[j]);
634 }
635 else
636 {
637 output->GetCellData()->AddArray(this->CellDataArray[j]);
638 }
639 double tempRange[2];
640 this->CellDataArray[j]->GetRange(tempRange, -1);
641 this->Minimum->InsertValue( j, (float)tempRange[0]);
642 this->Maximum->InsertValue( j, (float)tempRange[1]);
643 this->VectorLength->InsertValue( j, 1);
644 first = 1;
645 }
646 }
647 }
648
649
650 //----------------------------------------------------------------------------
RequestInformation(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * outputVector)651 int vtkMFIXReader::RequestInformation(
652 vtkInformation *vtkNotUsed(request),
653 vtkInformationVector **vtkNotUsed(inputVector),
654 vtkInformationVector *outputVector)
655 {
656 if ( this->RequestInformationFlag == 0)
657 {
658 if ( !this->FileName )
659 {
660 this->NumberOfPoints = 0;
661 this->NumberOfCells = 0;
662 vtkErrorMacro("No filename specified");
663 return 0;
664 }
665
666 this->SetProjectName(this->FileName);
667 this->ReadRestartFile();
668 this->CreateVariableNames();
669 this->GetTimeSteps();
670 this->CalculateMaxTimeStep();
671 this->MakeTimeStepTable(VariableNames->GetMaxId()+1);
672 this->GetNumberOfVariablesInSPXFiles();
673 this->MakeSPXTimeStepIndexTable(VariableNames->GetMaxId()+1);
674
675 for (int j = 0; j <= this->VariableNames->GetMaxId(); j++)
676 {
677 this->CellDataArraySelection->AddArray(
678 this->VariableNames->GetValue(j));
679 }
680
681 this->NumberOfPoints = (this->IMaximum2+1)
682 *(this->JMaximum2+1)*(this->KMaximum2+1);
683 this->NumberOfCells = this->IJKMaximum2;
684 this->NumberOfCellFields = this->VariableNames->GetMaxId()+1;
685 this->NumberOfTimeSteps = this->MaximumTimestep;
686 this->TimeStepRange[0] = 0;
687 this->TimeStepRange[1] = this->NumberOfTimeSteps-1;
688 this->RequestInformationFlag = 1;
689 this->GetAllTimes(outputVector);
690 }
691 return 1;
692 }
693
694 //----------------------------------------------------------------------------
GetNumberOfCellArrays()695 int vtkMFIXReader::GetNumberOfCellArrays()
696 {
697 return this->CellDataArraySelection->GetNumberOfArrays();
698 }
699
700 //----------------------------------------------------------------------------
GetCellArrayName(int index)701 const char* vtkMFIXReader::GetCellArrayName(int index)
702 {
703 return this->CellDataArraySelection->GetArrayName(index);
704 }
705
706 //----------------------------------------------------------------------------
GetCellArrayStatus(const char * name)707 int vtkMFIXReader::GetCellArrayStatus(const char* name)
708 {
709 return this->CellDataArraySelection->ArrayIsEnabled(name);
710 }
711
712 //----------------------------------------------------------------------------
SetCellArrayStatus(const char * name,int status)713 void vtkMFIXReader::SetCellArrayStatus(const char* name, int status)
714 {
715 if(status)
716 {
717 this->CellDataArraySelection->EnableArray(name);
718 }
719 else
720 {
721 this->CellDataArraySelection->DisableArray(name);
722 }
723 }
724
725 //----------------------------------------------------------------------------
DisableAllCellArrays()726 void vtkMFIXReader::DisableAllCellArrays()
727 {
728 this->CellDataArraySelection->DisableAllArrays();
729 }
730
731 //----------------------------------------------------------------------------
EnableAllCellArrays()732 void vtkMFIXReader::EnableAllCellArrays()
733 {
734 this->CellDataArraySelection->EnableAllArrays();
735 }
736
737 #if !defined(VTK_LEGACY_REMOVE)
738 //----------------------------------------------------------------------------
GetCellDataRange(int cellComp,int,float * min,float * max)739 void vtkMFIXReader::GetCellDataRange(int cellComp, int /* index */,
740 float *min, float *max)
741 {
742 #if !defined(VTK_LEGACY_SILENT)
743 vtkGenericWarningMacro("vtkMFIXReader::GetCellDataRange with \"index\" was deprecated in VTK 6.0");
744 #endif
745
746 *min = this->Minimum->GetValue(cellComp);
747 *max = this->Maximum->GetValue(cellComp);
748 }
749 #endif
750
751 //----------------------------------------------------------------------------
GetCellDataRange(int cellComp,float * min,float * max)752 void vtkMFIXReader::GetCellDataRange(int cellComp, float *min, float *max)
753 {
754 *min = this->Minimum->GetValue(cellComp);
755 *max = this->Maximum->GetValue(cellComp);
756 }
757
758 //----------------------------------------------------------------------------
SetProjectName(const char * infile)759 void vtkMFIXReader::SetProjectName (const char *infile)
760 {
761 int len = static_cast<int>(strlen(infile));
762 strncpy(this->RunName, infile, len-4);
763 }
764
765 //----------------------------------------------------------------------------
RestartVersionNumber(const char * buffer)766 void vtkMFIXReader::RestartVersionNumber(const char* buffer)
767 {
768 char s1[512];
769 char s2[512];
770 sscanf(buffer,"%s %s %f", s1, s2, &this->VersionNumber);
771 strncpy(this->Version, buffer, 100);
772 }
773
774 //----------------------------------------------------------------------------
GetInt(istream & in,int & val)775 void vtkMFIXReader::GetInt(istream& in, int &val)
776 {
777 in.read( (char*)&val,sizeof(int));
778 this->SwapInt(val);
779 }
780
781 //----------------------------------------------------------------------------
SwapInt(int & value)782 void vtkMFIXReader::SwapInt(int &value)
783 {
784 int result = ((value & 0x00FF) << 24) |
785 ((value & 0xFF00) << 8) |
786 ((value >> 8) & 0xFF00) |
787 ((value >> 24) & 0x00FF);
788 value = result;
789 }
790
791 //----------------------------------------------------------------------------
SwapDouble(double & value)792 void vtkMFIXReader::SwapDouble(double &value)
793 {
794 union Swap
795 {
796 double valDouble;
797 unsigned char valByte[8];
798 };
799
800 Swap source;
801 source.valDouble = value;
802
803 Swap result;
804 result.valByte[0] = source.valByte[7];
805 result.valByte[1] = source.valByte[6];
806 result.valByte[2] = source.valByte[5];
807 result.valByte[3] = source.valByte[4];
808 result.valByte[4] = source.valByte[3];
809 result.valByte[5] = source.valByte[2];
810 result.valByte[6] = source.valByte[1];
811 result.valByte[7] = source.valByte[0];
812
813 value = result.valDouble;
814 }
815
816 //----------------------------------------------------------------------------
SwapFloat(float & value)817 void vtkMFIXReader::SwapFloat(float &value)
818 {
819 union Swap
820 {
821 float valFloat;
822 unsigned char valByte[4];
823 };
824
825 Swap source;
826 source.valFloat = value;
827
828 Swap result;
829 result.valByte[0] = source.valByte[3];
830 result.valByte[1] = source.valByte[2];
831 result.valByte[2] = source.valByte[1];
832 result.valByte[3] = source.valByte[0];
833
834 value = result.valFloat;
835 }
836
837 //----------------------------------------------------------------------------
GetDouble(istream & in,double & val)838 void vtkMFIXReader::GetDouble(istream& in, double& val)
839 {
840 in.read( (char*)&val,sizeof(double));
841 this->SwapDouble(val);
842 }
843
844 //----------------------------------------------------------------------------
SkipBytes(istream & in,int n)845 void vtkMFIXReader::SkipBytes(istream& in, int n)
846 {
847 in.read(this->DataBuffer,n); // maybe seekg instead
848 }
849
850 //----------------------------------------------------------------------------
GetBlockOfDoubles(istream & in,vtkDoubleArray * v,int n)851 void vtkMFIXReader::GetBlockOfDoubles(istream& in, vtkDoubleArray *v, int n)
852 {
853 const int numberOfDoublesInBlock = 512/sizeof(double);
854 double tempArray[numberOfDoublesInBlock];
855 int numberOfRecords;
856
857 if ( n%numberOfDoublesInBlock == 0)
858 {
859 numberOfRecords = n/numberOfDoublesInBlock;
860 }
861 else
862 {
863 numberOfRecords = 1 + n/numberOfDoublesInBlock;
864 }
865
866 int c = 0;
867 for (int i=0; i<numberOfRecords; ++i)
868 {
869 in.read( (char*)&tempArray , 512 );
870 for (int j=0; j<numberOfDoublesInBlock; ++j)
871 {
872 if (c < n)
873 {
874 double temp = tempArray[j];
875 this->SwapDouble(temp);
876 v->InsertValue( c, temp);
877 ++c;
878 }
879 }
880 }
881 }
882
883 //----------------------------------------------------------------------------
GetBlockOfInts(istream & in,vtkIntArray * v,int n)884 void vtkMFIXReader::GetBlockOfInts(istream& in, vtkIntArray *v, int n)
885 {
886 const int numberOfIntsInBlock = 512/sizeof(int);
887 int tempArray[numberOfIntsInBlock];
888 int numberOfRecords;
889
890 if ( n%numberOfIntsInBlock == 0)
891 {
892 numberOfRecords = n/numberOfIntsInBlock;
893 }
894 else
895 {
896 numberOfRecords = 1 + n/numberOfIntsInBlock;
897 }
898
899 int c = 0;
900 for (int i = 0; i < numberOfRecords; ++i)
901 {
902 in.read( (char*)&tempArray , 512 );
903 for (int j=0; j<numberOfIntsInBlock; ++j)
904 {
905 if (c < n)
906 {
907 int temp = tempArray[j];
908 this->SwapInt(temp);
909 v->InsertValue( c, temp);
910 ++c;
911 }
912 }
913 }
914 }
915
916 //----------------------------------------------------------------------------
GetBlockOfFloats(istream & in,vtkFloatArray * v,int n)917 void vtkMFIXReader::GetBlockOfFloats(istream& in, vtkFloatArray *v, int n)
918 {
919 const int numberOfFloatsInBlock = 512/sizeof(float);
920 float tempArray[numberOfFloatsInBlock];
921 int numberOfRecords;
922
923 if ( n%numberOfFloatsInBlock == 0)
924 {
925 numberOfRecords = n/numberOfFloatsInBlock;
926 }
927 else
928 {
929 numberOfRecords = 1 + n/numberOfFloatsInBlock;
930 }
931
932 bool modified = false;
933 int c = 0;
934 int cnt = 0;
935 for (int i=0; i<numberOfRecords; ++i)
936 {
937 in.read( (char*)&tempArray , 512 );
938 for (int j=0; j<numberOfFloatsInBlock; ++j)
939 {
940 if (c < n)
941 {
942 float temp = tempArray[j];
943 this->SwapFloat(temp);
944 if ( this->Flag->GetValue(c) < 10)
945 {
946 v->InsertValue(cnt, temp);
947 cnt++;
948 modified = true;
949 }
950 ++c;
951 }
952 }
953 }
954
955 if (modified)
956 {
957 v->Modified();
958 }
959 }
960
961 //----------------------------------------------------------------------------
ReadRestartFile()962 void vtkMFIXReader::ReadRestartFile()
963 {
964 int dimensionUsr = 5;
965
966 #ifdef _WIN32
967 ifstream in(this->FileName,ios::binary);
968 #else
969 ifstream in(this->FileName);
970 #endif
971
972 if (!in)
973 {
974 //cout << "could not open file" << endl;
975 return;
976 }
977
978 this->DataBuffer[512] = '\0';
979
980 // version : record 1
981 memset(this->DataBuffer,0,513);
982 in.read(this->DataBuffer,512);
983 this->RestartVersionNumber(this->DataBuffer);
984
985 // skip 2 linesline : records 2 and 3
986 in.read(this->DataBuffer,512);
987 in.read(this->DataBuffer,512);
988
989 // IMinimum1 etc : record 4
990 memset(this->DataBuffer,0,513);
991
992 if (strcmp(this->Version, "RES = 01.00") == 0)
993 {
994 this->GetInt(in,this->IMinimum1);
995 this->GetInt(in,this->JMinimum1);
996 this->GetInt(in,this->KMinimum1);
997 this->GetInt(in,this->IMaximum);
998 this->GetInt(in,this->JMaximum);
999 this->GetInt(in,this->KMaximum);
1000 this->GetInt(in,this->IMaximum1);
1001 this->GetInt(in,this->JMaximum1);
1002 this->GetInt(in,this->KMaximum1);
1003 this->GetInt(in,this->IMaximum2);
1004 this->GetInt(in,this->JMaximum2);
1005 this->GetInt(in,this->KMaximum2);
1006 this->GetInt(in,this->IJMaximum2);
1007 this->GetInt(in,this->IJKMaximum2);
1008 this->GetInt(in,this->MMAX);
1009 this->GetDouble(in,this->DeltaTime);
1010 this->GetDouble(in,this->XLength);
1011 this->GetDouble(in,this->YLength);
1012 this->GetDouble(in,this->ZLength);
1013
1014 // 15 ints ... 4 doubles = 92 bytes
1015 this->SkipBytes(in,420);
1016 }
1017 else if (strcmp(this->Version, "RES = 01.01") == 0 ||
1018 strcmp(this->Version, "RES = 01.02") == 0)
1019 {
1020 this->GetInt(in,this->IMinimum1);
1021 this->GetInt(in,this->JMinimum1);
1022 this->GetInt(in,this->KMinimum1);
1023 this->GetInt(in,this->IMaximum);
1024 this->GetInt(in,this->JMaximum);
1025 this->GetInt(in,this->KMaximum);
1026 this->GetInt(in,this->IMaximum1);
1027 this->GetInt(in,this->JMaximum1);
1028 this->GetInt(in,this->KMaximum1);
1029 this->GetInt(in,this->IMaximum2);
1030 this->GetInt(in,this->JMaximum2);
1031 this->GetInt(in,this->KMaximum2);
1032 this->GetInt(in,this->IJMaximum2);
1033 this->GetInt(in,this->IJKMaximum2);
1034 this->GetInt(in,this->MMAX);
1035 this->GetInt(in,this->DimensionIc);
1036 this->GetInt(in,this->DimensionBc);
1037 this->GetDouble(in,this->DeltaTime);
1038 this->GetDouble(in,this->XLength);
1039 this->GetDouble(in,this->YLength);
1040 this->GetDouble(in,this->ZLength);
1041
1042 // 17 ints ... 4 doubles = 100 bytes
1043 this->SkipBytes(in,412);
1044 }
1045 else if(strcmp(this->Version, "RES = 01.03") == 0)
1046 {
1047 this->GetInt(in,this->IMinimum1);
1048 this->GetInt(in,this->JMinimum1);
1049 this->GetInt(in,this->KMinimum1);
1050 this->GetInt(in,this->IMaximum);
1051 this->GetInt(in,this->JMaximum);
1052 this->GetInt(in,this->KMaximum);
1053 this->GetInt(in,this->IMaximum1);
1054 this->GetInt(in,this->JMaximum1);
1055 this->GetInt(in,this->KMaximum1);
1056 this->GetInt(in,this->IMaximum2);
1057 this->GetInt(in,this->JMaximum2);
1058 this->GetInt(in,this->KMaximum2);
1059 this->GetInt(in,this->IJMaximum2);
1060 this->GetInt(in,this->IJKMaximum2);
1061 this->GetInt(in,this->MMAX);
1062 this->GetInt(in,this->DimensionIc);
1063 this->GetInt(in,this->DimensionBc);
1064 this->GetDouble(in,this->DeltaTime);
1065 this->GetDouble(in,this->XMinimum);
1066 this->GetDouble(in,this->XLength);
1067 this->GetDouble(in,this->YLength);
1068 this->GetDouble(in,this->ZLength);
1069
1070 // 17 ints ... 5 doubles = 108 bytes
1071 this->SkipBytes(in,404);
1072 }
1073 else if(strcmp(this->Version, "RES = 01.04") == 0)
1074 {
1075 this->GetInt(in,this->IMinimum1);
1076 this->GetInt(in,this->JMinimum1);
1077 this->GetInt(in,this->KMinimum1);
1078 this->GetInt(in,this->IMaximum);
1079 this->GetInt(in,this->JMaximum);
1080 this->GetInt(in,this->KMaximum);
1081 this->GetInt(in,this->IMaximum1);
1082 this->GetInt(in,this->JMaximum1);
1083 this->GetInt(in,this->KMaximum1);
1084 this->GetInt(in,this->IMaximum2);
1085 this->GetInt(in,this->JMaximum2);
1086 this->GetInt(in,this->KMaximum2);
1087 this->GetInt(in,this->IJMaximum2);
1088 this->GetInt(in,this->IJKMaximum2);
1089 this->GetInt(in,this->MMAX);
1090 this->GetInt(in,this->DimensionIc);
1091 this->GetInt(in,this->DimensionBc);
1092 this->GetInt(in,this->DimensionC);
1093 this->GetDouble(in,this->DeltaTime);
1094 this->GetDouble(in,this->XMinimum);
1095 this->GetDouble(in,this->XLength);
1096 this->GetDouble(in,this->YLength);
1097 this->GetDouble(in,this->ZLength);
1098
1099 // 18 ints ... 5 doubles = 112 bytes
1100 this->SkipBytes(in,400);
1101 }
1102 else if(strcmp(this->Version, "RES = 01.05") == 0)
1103 {
1104 this->GetInt(in,this->IMinimum1);
1105 this->GetInt(in,this->JMinimum1);
1106 this->GetInt(in,this->KMinimum1);
1107 this->GetInt(in,this->IMaximum);
1108 this->GetInt(in,this->JMaximum);
1109 this->GetInt(in,this->KMaximum);
1110 this->GetInt(in,this->IMaximum1);
1111 this->GetInt(in,this->JMaximum1);
1112 this->GetInt(in,this->KMaximum1);
1113 this->GetInt(in,this->IMaximum2);
1114 this->GetInt(in,this->JMaximum2);
1115 this->GetInt(in,this->KMaximum2);
1116 this->GetInt(in,this->IJMaximum2);
1117 this->GetInt(in,this->IJKMaximum2);
1118 this->GetInt(in,this->MMAX);
1119 this->GetInt(in,this->DimensionIc);
1120 this->GetInt(in,this->DimensionBc);
1121 this->GetInt(in,this->DimensionC);
1122 this->GetInt(in,this->DimensionIs);
1123 this->GetDouble(in,this->DeltaTime);
1124 this->GetDouble(in,this->XMinimum);
1125 this->GetDouble(in,this->XLength);
1126 this->GetDouble(in,this->YLength);
1127 this->GetDouble(in,this->ZLength);
1128
1129 // 19 ints ... 5 doubles = 116 bytes
1130 this->SkipBytes(in,396);
1131 }
1132 else
1133 {
1134 this->GetInt(in,this->IMinimum1);
1135 this->GetInt(in,this->JMinimum1);
1136 this->GetInt(in,this->KMinimum1);
1137 this->GetInt(in,this->IMaximum);
1138 this->GetInt(in,this->JMaximum);
1139 this->GetInt(in,this->KMaximum);
1140 this->GetInt(in,this->IMaximum1);
1141 this->GetInt(in,this->JMaximum1);
1142 this->GetInt(in,this->KMaximum1);
1143 this->GetInt(in,this->IMaximum2);
1144 this->GetInt(in,this->JMaximum2);
1145 this->GetInt(in,this->KMaximum2);
1146 this->GetInt(in,this->IJMaximum2);
1147 this->GetInt(in,this->IJKMaximum2);
1148 this->GetInt(in,this->MMAX);
1149 this->GetInt(in,this->DimensionIc);
1150 this->GetInt(in,this->DimensionBc);
1151 this->GetInt(in,this->DimensionC);
1152 this->GetInt(in,this->DimensionIs);
1153 this->GetDouble(in,this->DeltaTime);
1154 this->GetDouble(in,this->XMinimum);
1155 this->GetDouble(in,this->XLength);
1156 this->GetDouble(in,this->YLength);
1157 this->GetDouble(in,this->ZLength);
1158 this->GetDouble(in,this->Ce);
1159 this->GetDouble(in,this->Cf);
1160 this->GetDouble(in,this->Phi);
1161 this->GetDouble(in,this->PhiW);
1162
1163 // 19 ints ... 9 doubles = 148 bytes
1164 this->SkipBytes(in,364);
1165 }
1166
1167 const int numberOfFloatsInBlock = 512/sizeof(float);
1168
1169 if ( this->IJKMaximum2%numberOfFloatsInBlock == 0)
1170 {
1171 this->SPXRecordsPerTimestep = this->IJKMaximum2/numberOfFloatsInBlock;
1172 }
1173 else
1174 {
1175 this->SPXRecordsPerTimestep = 1 + this->IJKMaximum2/numberOfFloatsInBlock;
1176 }
1177
1178 // C , C_name and nmax
1179
1180 this->NMax->Resize(this->MMAX+1);
1181 for (int lc=0; lc<this->MMAX+1; ++lc)
1182 {
1183 this->NMax->InsertValue(lc, 1);
1184 }
1185
1186 this->C->Resize(this->DimensionC);
1187
1188 if (this->VersionNumber > 1.04)
1189 {
1190 this->GetBlockOfDoubles (in, this->C, this->DimensionC);
1191
1192 for (int lc=0; lc<DimensionC; ++lc)
1193 {
1194 in.read(this->DataBuffer,512); // c_name[]
1195 }
1196
1197 if (this->VersionNumber < 1.12)
1198 {
1199 this->GetBlockOfInts(in, this->NMax,this->MMAX+1);
1200 }
1201 else
1202 {
1203 // what is the diff between this and above ???
1204 for (int lc=0; lc<this->MMAX+1; ++lc)
1205 {
1206 int temp;
1207 this->GetInt(in,temp);
1208 this->NMax->InsertValue(lc, temp);
1209 }
1210 this->SkipBytes(in,512-(this->MMAX+1)*sizeof(int));
1211 }
1212 }
1213
1214 this->Dx->Resize(this->IMaximum2);
1215 this->Dy->Resize(this->JMaximum2);
1216 this->Dz->Resize(this->KMaximum2);
1217
1218 this->GetBlockOfDoubles(in, this->Dx,this->IMaximum2);
1219 this->GetBlockOfDoubles(in, this->Dy,this->JMaximum2);
1220 this->GetBlockOfDoubles(in, this->Dz,this->KMaximum2);
1221
1222 // RunName etc.
1223
1224 memset(this->Units,0,17);
1225 memset(this->CoordinateSystem,0,17);
1226
1227 in.read(this->DataBuffer,120); // run_name , description
1228 in.read(this->Units,16); // Units
1229 in.read(this->DataBuffer,16); // run_type
1230 in.read(this->CoordinateSystem,16); // CoordinateSystem
1231
1232 this->SkipBytes(in,512-168);
1233
1234 char tempCharArray[17];
1235
1236 memset(tempCharArray,0,17);
1237
1238 int ic = 0;
1239 for (int i=0; i<17; ++i)
1240 {
1241 if (this->Units[i] != ' ')
1242 {
1243 tempCharArray[ic++] = this->Units[i];
1244 }
1245 }
1246
1247 memset(tempCharArray,0,17);
1248
1249 ic = 0;
1250 for (int i=0; i<17; ++i)
1251 {
1252 if (this->CoordinateSystem[i] != ' ')
1253 {
1254 tempCharArray[ic++] = this->CoordinateSystem[i];
1255 }
1256 }
1257 strcpy(this->CoordinateSystem,tempCharArray);
1258
1259 if (this->VersionNumber >= 1.04)
1260 {
1261 this->TempD->Resize(this->NMax->GetValue(0));
1262 this->GetBlockOfDoubles(in, this->TempD, this->NMax->GetValue(0)); // MW_g
1263 for (int i=0; i<this->MMAX; ++i)
1264 {
1265 in.read(this->DataBuffer,512); // MW_s
1266 }
1267 }
1268
1269 in.read(this->DataBuffer,512); // D_p etc.
1270
1271 // read in the "DimensionIc" variables (and ignore ... not used by ani_mfix)
1272 this->TempI->Resize(this->DimensionIc);
1273 this->TempD->Resize(this->DimensionIc);
1274
1275 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_x_w
1276 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_x_e
1277 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_y_s
1278 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_y_n
1279 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_z_b
1280 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_z_t
1281
1282 this->GetBlockOfInts(in, this->TempI,this->DimensionIc); // ic_i_w
1283 this->GetBlockOfInts(in, this->TempI,this->DimensionIc); // ic_i_e
1284 this->GetBlockOfInts(in, this->TempI,this->DimensionIc); // ic_j_s
1285 this->GetBlockOfInts(in, this->TempI,this->DimensionIc); // ic_j_n
1286 this->GetBlockOfInts(in, this->TempI,this->DimensionIc); // ic_k_b
1287 this->GetBlockOfInts(in, this->TempI,this->DimensionIc); // ic_k_t
1288
1289 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_ep_g
1290 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_p_g
1291 this->GetBlockOfDoubles(in, this->TempD,this->DimensionIc); // ic_t_g
1292
1293 if (this->VersionNumber < 1.15)
1294 {
1295 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_t_s(1,1)
1296 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_t_s(1,2)
1297 // or ic_tmp
1298 }
1299
1300 if (this->VersionNumber >= 1.04)
1301 {
1302 for (int i=0; i<this->NMax->GetValue(0); ++i)
1303 {
1304 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_x_g
1305 }
1306 }
1307
1308 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_u_g
1309 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_v_g
1310 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_w_g
1311
1312 for (int lc=0; lc<this->MMAX; ++lc)
1313 {
1314 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_rop_s
1315 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_u_s
1316 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_v_s
1317 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_w_s
1318
1319 if (this->VersionNumber >= 1.15)
1320 {
1321 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_t_s
1322 }
1323
1324 if (this->VersionNumber >= 1.04)
1325 {
1326 for (int n=0; n<this->NMax->GetValue(lc+1); ++n)
1327 {
1328 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_x_s
1329 }
1330 }
1331 }
1332
1333 // read in the "DimensionBc" variables (and ignore ... not used by ani_mfix)
1334 this->TempI->Resize(this->DimensionBc);
1335 this->TempD->Resize(this->DimensionBc);
1336
1337 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_x_w
1338 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_x_e
1339 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc y s
1340 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc y n
1341 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc z b
1342 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc z t
1343 this->GetBlockOfInts(in,this->TempI,this->DimensionBc); // bc i w
1344 this->GetBlockOfInts(in,this->TempI,this->DimensionBc); // bc i e
1345 this->GetBlockOfInts(in,this->TempI,this->DimensionBc); // bc j s
1346 this->GetBlockOfInts(in,this->TempI,this->DimensionBc); // bc j n
1347 this->GetBlockOfInts(in,this->TempI,this->DimensionBc); // bc k b
1348 this->GetBlockOfInts(in,this->TempI,this->DimensionBc); // bc k t
1349 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc ep g
1350 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc p g
1351 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc t g
1352
1353 if (VersionNumber < 1.15)
1354 {
1355 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_t_s(1,1)
1356 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_t_s(1,1)
1357 // or bc_tmp
1358 }
1359
1360 if (VersionNumber >= 1.04)
1361 {
1362 for (int i=0; i<this->NMax->GetValue(0); ++i)
1363 {
1364 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_x_g
1365 }
1366 }
1367
1368 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc u g
1369 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc v g
1370 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc w g
1371 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc ro g
1372 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_rop_g
1373 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc volflow g
1374 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc massflow g
1375
1376 for (int lc=0; lc<this->MMAX; ++lc)
1377 {
1378 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc rop s
1379 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc u s
1380 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc v s
1381
1382 if (this->VersionNumber >= 1.04)
1383 {
1384 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc w s
1385
1386 if (this->VersionNumber >= 1.15)
1387 {
1388 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc t s
1389 }
1390 for (int n=0; n<this->NMax->GetValue(lc+1); ++n)
1391 {
1392 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc x s
1393 }
1394 }
1395 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc volflow s
1396 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc massflow s
1397 }
1398
1399 if (strcmp(this->Version,"RES = 01.00") == 0)
1400 {
1401 for (int lc=0; lc<10; ++lc)
1402 {
1403 in.read(this->DataBuffer,512); // BC TYPE
1404 }
1405 }
1406 else
1407 {
1408 for (int lc=0; lc<this->DimensionBc; ++lc)
1409 {
1410 in.read(this->DataBuffer,512); // BC TYPE
1411 }
1412 }
1413
1414 this->Flag->Resize(this->IJKMaximum2);
1415 this->GetBlockOfInts(in, this->Flag,this->IJKMaximum2);
1416
1417 // DimensionIs varibles (not needed by ani_mfix)
1418 this->TempI->Resize(this->DimensionIs);
1419 this->TempD->Resize(this->DimensionIs);
1420
1421 if (this->VersionNumber >= 1.04)
1422 {
1423 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is x w
1424 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is x e
1425 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is y s
1426 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is y n
1427 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is z b
1428 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is z t
1429 this->GetBlockOfInts(in,this->TempI,this->DimensionIs); // is i w
1430 this->GetBlockOfInts(in,this->TempI,this->DimensionIs); // is i e
1431 this->GetBlockOfInts(in,this->TempI,this->DimensionIs); // is j s
1432 this->GetBlockOfInts(in,this->TempI,this->DimensionIs); // is j n
1433 this->GetBlockOfInts(in,this->TempI,this->DimensionIs); // is k b
1434 this->GetBlockOfInts(in,this->TempI,this->DimensionIs); // is k t
1435 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is_pc(1,1)
1436 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs); // is_pc(1,2)
1437
1438 if (this->VersionNumber >= 1.07)
1439 {
1440 for (int i=0; i<this->MMAX; ++i)
1441 {
1442 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIs);//is_vel_s
1443 }
1444 }
1445
1446 for (int lc=0; lc<this->DimensionIs; ++lc)
1447 {
1448 in.read(this->DataBuffer,512); // is_type
1449 }
1450 }
1451
1452 if (this->VersionNumber >= 1.08)
1453 {
1454 in.read(this->DataBuffer,512);
1455 }
1456
1457 if (this->VersionNumber >= 1.09)
1458 {
1459 in.read(this->DataBuffer,512);
1460
1461 if (this->VersionNumber >= 1.5)
1462 {
1463 this->GetInt(in,this->NumberOfSPXFilesUsed);
1464 this->SkipBytes(in,508);
1465 }
1466
1467 for (int lc=0; lc< this->NumberOfSPXFilesUsed; ++lc)
1468 {
1469 in.read(this->DataBuffer,512); // spx_dt
1470 }
1471
1472 for (int lc=0; lc<this->MMAX+1; ++lc)
1473 {
1474 in.read(this->DataBuffer,512); // species_eq
1475 }
1476
1477 this->TempD->Resize(dimensionUsr);
1478
1479 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr_dt
1480 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr x w
1481 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr x e
1482 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr y s
1483 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr y n
1484 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr z b
1485 this->GetBlockOfDoubles(in,this->TempD,dimensionUsr); // usr z t
1486
1487 for (int lc=0; lc<dimensionUsr; ++lc)
1488 {
1489 in.read(this->DataBuffer,512); // usr_ext etc.
1490 }
1491
1492 this->TempD->Resize(this->DimensionIc);
1493 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_p_star
1494 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_l_scale
1495 for (int lc=0; lc<this->DimensionIc; ++lc)
1496 {
1497 in.read(this->DataBuffer,512); // ic_type
1498 }
1499
1500 this->TempD->Resize(DimensionBc);
1501 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_dt_0
1502 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_jet_g0
1503 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_dt_h
1504 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_jet_gh
1505 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_dt_l
1506 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_jet_gl
1507 }
1508
1509 if (this->VersionNumber >= 1.1)
1510 {
1511 in.read(this->DataBuffer,512); // mu_gmax
1512 }
1513
1514 if (this->VersionNumber >= 1.11)
1515 {
1516 in.read(this->DataBuffer,512); // x_ex , model_b
1517 }
1518
1519 if (this->VersionNumber >= 1.12)
1520 {
1521 in.read(this->DataBuffer,512); // p_ref , etc.
1522 in.read(this->DataBuffer,512); // leq_it , leq_method
1523
1524 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_hw_g
1525 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_uw_g
1526 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_vw_g
1527 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_ww_g
1528
1529 for (int lc=0; lc<this->MMAX; ++lc)
1530 {
1531 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_hw_s
1532 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_uw_s
1533 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_vw_s
1534 this->GetBlockOfDoubles(in,this->TempD,this->DimensionBc); // bc_ww_s
1535 }
1536 }
1537
1538 if (this->VersionNumber >= 1.13)
1539 {
1540 in.read(this->DataBuffer,512); // momentum_x_eq , etc.
1541 }
1542
1543 if (this->VersionNumber >= 1.14)
1544 {
1545 in.read(this->DataBuffer,512); // detect_small
1546 }
1547
1548 if (this->VersionNumber >= 1.15)
1549 {
1550 in.read(this->DataBuffer,512); // k_g0 , etc.
1551
1552 this->TempD->Resize(this->DimensionIc);
1553
1554 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_gama_rg
1555 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_t_rg
1556
1557 for (int lc=0; lc<this->MMAX; ++lc)
1558 {
1559 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_gama_rs
1560 this->GetBlockOfDoubles(in,this->TempD,this->DimensionIc); // ic_t_rs
1561 }
1562 }
1563
1564 if (this->VersionNumber >= 1.2)
1565 {
1566 in.read(this->DataBuffer,512); // norm_g , norm_s
1567 }
1568
1569 if (this->VersionNumber >= 1.3)
1570 {
1571 this->GetInt(in,this->NumberOfScalars);
1572 this->SkipBytes(in,sizeof(double)); // tol_resid_scalar
1573
1574 int DIM_tmp;
1575 this->GetInt(in,DIM_tmp);
1576 this->SkipBytes(in,512-sizeof(double)-2*sizeof(int));
1577
1578 this->TempI->Resize(DIM_tmp);
1579 this->GetBlockOfInts(in,this->TempI,DIM_tmp); // Phase4Scalar;
1580 }
1581
1582 if (this->VersionNumber >= 1.5)
1583 {
1584 this->GetInt(in,this->NumberOfReactionRates);
1585 this->SkipBytes(in,508);
1586 }
1587
1588 if (this->VersionNumber >= 1.5999)
1589 {
1590 int tmp;
1591 this->GetInt(in,tmp);
1592 this->SkipBytes(in,508);
1593
1594 if (tmp != 0)
1595 {
1596 this->BkEpsilon = true;
1597 }
1598 }
1599 if (this->VersionNumber >= 1.7999)
1600 {
1601 for( int i = 0; i < this->MMAX; ++i)
1602 {
1603 this->SkipBytes(in,512);
1604 }
1605 }
1606 in.close();
1607 }
1608
1609 //----------------------------------------------------------------------------
CreateVariableNames()1610 void vtkMFIXReader::CreateVariableNames()
1611 {
1612 char fileName[256];
1613 int cnt = 0;
1614 char uString[120];
1615 char vString[120];
1616 char wString[120];
1617 char svString[120];
1618 char tempString[120];
1619 char ropString[120];
1620 char temperatureString[120];
1621 char variableString[120];
1622
1623 for (int i=0; i<this->NumberOfSPXFilesUsed; ++i)
1624 {
1625 for(int k = 0; k < (int)sizeof(fileName); k++)
1626 {
1627 fileName[k]=0;
1628 }
1629 strncpy(fileName, this->FileName, strlen(this->FileName)-4);
1630
1631 if (i==0)
1632 {
1633 strcat(fileName, ".SP1");
1634 }
1635 else if (i==1)
1636 {
1637 strcat(fileName, ".SP2");
1638 }
1639 else if (i==2)
1640 {
1641 strcat(fileName, ".SP3");
1642 }
1643 else if (i==3)
1644 {
1645 strcat(fileName, ".SP4");
1646 }
1647 else if (i==4)
1648 {
1649 strcat(fileName, ".SP5");
1650 }
1651 else if (i==5)
1652 {
1653 strcat(fileName, ".SP6");
1654 }
1655 else if (i==6)
1656 {
1657 strcat(fileName, ".SP7");
1658 }
1659 else if (i==7)
1660 {
1661 strcat(fileName, ".SP8");
1662 }
1663 else if (i==8)
1664 {
1665 strcat(fileName, ".SP9");
1666 }
1667 else if (i==9)
1668 {
1669 strcat(fileName, ".SPA");
1670 }
1671 else
1672 {
1673 strcat(fileName, ".SPB");
1674 }
1675
1676 #ifdef _WIN32
1677 ifstream in(fileName,ios::binary);
1678 #else
1679 ifstream in(fileName);
1680 #endif
1681 if (in) // file exists
1682 {
1683 this->SpxFileExists->InsertValue(i, 1);
1684
1685 switch (i+1)
1686 {
1687
1688 case 1:
1689 this->VariableNames->InsertValue(cnt++,"EP_g");
1690 this->VariableIndexToSPX->InsertValue(cnt-1, 1);
1691 this->VariableComponents->InsertValue(cnt-1, 1);
1692 break;
1693
1694 case 2:
1695 this->VariableNames->InsertValue(cnt++,"P_g");
1696 this->VariableIndexToSPX->InsertValue(cnt-1, 2);
1697 this->VariableComponents->InsertValue(cnt-1, 1);
1698 this->VariableNames->InsertValue(cnt++,"P_star");
1699 this->VariableIndexToSPX->InsertValue(cnt-1, 2);
1700 this->VariableComponents->InsertValue(cnt-1, 1);
1701 break;
1702
1703 case 3:
1704 this->VariableNames->InsertValue(cnt++,"U_g");
1705 this->VariableIndexToSPX->InsertValue(cnt-1, 3);
1706 this->VariableComponents->InsertValue(cnt-1, 1);
1707 this->VariableNames->InsertValue(cnt++,"V_g");
1708 this->VariableIndexToSPX->InsertValue(cnt-1, 3);
1709 this->VariableComponents->InsertValue(cnt-1, 1);
1710 this->VariableNames->InsertValue(cnt++,"W_g");
1711 this->VariableIndexToSPX->InsertValue(cnt-1, 3);
1712 this->VariableComponents->InsertValue(cnt-1, 1);
1713 this->VariableNames->InsertValue(cnt++,"Gas Velocity");
1714 this->VariableIndexToSPX->InsertValue(cnt-1, 3);
1715 this->VariableComponents->InsertValue(cnt-1, 3);
1716 break;
1717
1718 case 4:
1719 for (int j=0; j<this->MMAX; ++j)
1720 {
1721 for(int k=0;k<(int)sizeof(uString);k++)
1722 {
1723 uString[k]=0;
1724 }
1725 for(int k=0;k<(int)sizeof(vString);k++)
1726 {
1727 vString[k]=0;
1728 }
1729 for(int k=0;k<(int)sizeof(wString);k++)
1730 {
1731 wString[k]=0;
1732 }
1733 for(int k=0;k<(int)sizeof(svString);k++)
1734 {
1735 svString[k]=0;
1736 }
1737 strcpy(uString, "U_s_");
1738 strcpy(vString, "V_s_");
1739 strcpy(wString, "W_s_");
1740 strcpy(svString, "Solids_Velocity_");
1741 sprintf(tempString, "%d", j+1);
1742 strcat(uString, tempString);
1743 strcat(vString, tempString);
1744 strcat(wString, tempString);
1745 strcat(svString, tempString);
1746 this->VariableNames->InsertValue(cnt++, uString);
1747 this->VariableIndexToSPX->InsertValue(cnt-1, 4);
1748 this->VariableComponents->InsertValue(cnt-1, 1);
1749
1750 this->VariableNames->InsertValue(cnt++, vString);
1751 this->VariableIndexToSPX->InsertValue(cnt-1, 4);
1752 this->VariableComponents->InsertValue(cnt-1, 1);
1753
1754 this->VariableNames->InsertValue(cnt++, wString);
1755 this->VariableIndexToSPX->InsertValue(cnt-1, 4);
1756 this->VariableComponents->InsertValue(cnt-1, 1);
1757
1758 this->VariableNames->InsertValue(cnt++, svString);
1759 this->VariableIndexToSPX->InsertValue(cnt-1, 4);
1760 this->VariableComponents->InsertValue(cnt-1, 3);
1761 }
1762 break;
1763
1764 case 5:
1765 for (int j=0; j<this->MMAX; ++j)
1766 {
1767 for(int k=0;k<(int)sizeof(ropString);k++)
1768 {
1769 ropString[k]=0;
1770 }
1771 strcpy(ropString, "ROP_s_");
1772 sprintf(tempString, "%d", j+1);
1773 strcat(ropString, tempString);
1774 this->VariableNames->InsertValue(cnt++, ropString);
1775 this->VariableIndexToSPX->InsertValue(cnt-1, 5);
1776 this->VariableComponents->InsertValue(cnt-1, 1);
1777 }
1778 break;
1779
1780 case 6:
1781 this->VariableNames->InsertValue(cnt++, "T_g");
1782 this->VariableIndexToSPX->InsertValue(cnt-1, 6);
1783 this->VariableComponents->InsertValue(cnt-1, 1);
1784
1785 if (this->VersionNumber <= 1.15)
1786 {
1787 this->VariableNames->InsertValue(cnt++, "T_s_1");
1788 this->VariableIndexToSPX->InsertValue(cnt-1, 6);
1789 this->VariableComponents->InsertValue(cnt-1, 1);
1790
1791 if (this->MMAX > 1)
1792 {
1793 this->VariableNames->InsertValue(cnt++, "T_s_2");
1794 this->VariableIndexToSPX->InsertValue(cnt-1, 6);
1795 this->VariableComponents->InsertValue(cnt-1, 1);
1796 }
1797 else
1798 {
1799 this->VariableNames->InsertValue(cnt++, "T_s_2_not_used");
1800 this->VariableIndexToSPX->InsertValue(cnt-1, 6);
1801 this->VariableComponents->InsertValue(cnt-1, 1);
1802 }
1803 }
1804 else
1805 {
1806 for (int j=0; j<this->MMAX; ++j)
1807 {
1808 for(int k=0;k<(int)sizeof(temperatureString);k++)
1809 {
1810 temperatureString[k]=0;
1811 }
1812 strcpy(temperatureString, "T_s_");
1813 sprintf(tempString, "%d", j+1);
1814 strcat(temperatureString, tempString);
1815 this->VariableNames->InsertValue(cnt++, temperatureString);
1816 this->VariableIndexToSPX->InsertValue(cnt-1, 6);
1817 this->VariableComponents->InsertValue(cnt-1, 1);
1818 }
1819 }
1820 break;
1821
1822 case 7:
1823 for (int j=0; j<this->NMax->GetValue(0); ++j)
1824 {
1825 for (int k=0;k<(int)sizeof(variableString);k++)
1826 {
1827 variableString[k]=0;
1828 }
1829 strcpy(variableString, "X_g_");
1830 sprintf(tempString, "%d", j+1);
1831 strcat(variableString, tempString);
1832 this->VariableNames->InsertValue(cnt++, variableString);
1833 this->VariableIndexToSPX->InsertValue(cnt-1, 7);
1834 this->VariableComponents->InsertValue(cnt-1, 1);
1835 }
1836
1837 for (int m=1; m<=this->MMAX; ++m)
1838 {
1839 for (int j=0; j<this->NMax->GetValue(m); ++j)
1840 {
1841 char tempString1[120];
1842 char tempString2[120];
1843 for (int k=0;k<(int)sizeof(variableString);k++)
1844 {
1845 variableString[k]=0;
1846 }
1847 strcpy(variableString, "X_s_");
1848 sprintf(tempString1, "%d", m);
1849 sprintf(tempString2, "%d", j+1);
1850 strcat(variableString, tempString1);
1851 strcat(variableString, "_");
1852 strcat(variableString, tempString2);
1853 this->VariableNames->InsertValue(cnt++, variableString);
1854 this->VariableIndexToSPX->InsertValue(cnt-1, 7);
1855 this->VariableComponents->InsertValue(cnt-1, 1);
1856 }
1857 }
1858 break;
1859
1860 case 8:
1861 for (int j=0; j<MMAX; ++j)
1862 {
1863 for (int k=0;k<(int)sizeof(variableString);k++)
1864 {
1865 variableString[k]=0;
1866 }
1867 strcpy(variableString, "Theta_m_");
1868 sprintf(tempString, "%d", j+1);
1869 strcat(variableString, tempString);
1870 this->VariableNames->InsertValue(cnt++, variableString);
1871 this->VariableIndexToSPX->InsertValue(cnt-1, 8);
1872 this->VariableComponents->InsertValue(cnt-1, 1);
1873 }
1874 break;
1875
1876 case 9:
1877 for (int j=0; j<this->NumberOfScalars; ++j)
1878 {
1879 for(int k=0;k<(int)sizeof(variableString);k++)
1880 {
1881 variableString[k]=0;
1882 }
1883 strcpy(variableString, "Scalar_");
1884 sprintf(tempString, "%d", j+1);
1885 strcat(variableString, tempString);
1886 this->VariableNames->InsertValue(cnt++, variableString);
1887 this->VariableIndexToSPX->InsertValue(cnt-1, 9);
1888 this->VariableComponents->InsertValue(cnt-1, 1);
1889 }
1890 break;
1891
1892 case 10:
1893 for (int j=0; j<this->NumberOfReactionRates; ++j)
1894 {
1895 for (int k=0;k<(int)sizeof(variableString);k++)
1896 {
1897 variableString[k]=0;
1898 }
1899 strcpy(variableString, "RRates_");
1900 sprintf(tempString, "%d", j+1);
1901 strcat(variableString, tempString);
1902 this->VariableNames->InsertValue(cnt++, variableString);
1903 this->VariableIndexToSPX->InsertValue(cnt-1, 10);
1904 this->VariableComponents->InsertValue(cnt-1, 1);
1905 }
1906 break;
1907
1908 case 11:
1909 if (BkEpsilon)
1910 {
1911 this->VariableNames->InsertValue(cnt++, "k_turb_g");
1912 this->VariableIndexToSPX->InsertValue(cnt-1, 11);
1913 this->VariableComponents->InsertValue(cnt-1, 1);
1914 this->VariableNames->InsertValue(cnt++, "e_turb_g");
1915 this->VariableIndexToSPX->InsertValue(cnt-1, 11);
1916 this->VariableComponents->InsertValue(cnt-1, 1);
1917 }
1918 break;
1919 default:
1920 vtkWarningMacro(<< "unknown SPx file : " << i << "\n");
1921 break;
1922 }
1923 }
1924 else
1925 {
1926 this->SpxFileExists->InsertValue(i, 0);
1927 }
1928 in.close();
1929 }
1930 }
1931
1932 //----------------------------------------------------------------------------
GetTimeSteps()1933 void vtkMFIXReader::GetTimeSteps()
1934 {
1935 int nextRecord, numberOfRecords;
1936 char fileName[256];
1937 int cnt = 0;
1938
1939 for (int i=0; i<this->NumberOfSPXFilesUsed; ++i)
1940 {
1941 for (int k=0;k<(int)sizeof(fileName);k++)
1942 {
1943 fileName[k]=0;
1944 }
1945 strncpy(fileName, this->FileName, strlen(this->FileName)-4);
1946 if (i==0)
1947 {
1948 strcat(fileName, ".SP1");
1949 }
1950 else if (i==1)
1951 {
1952 strcat(fileName, ".SP2");
1953 }
1954 else if (i==2)
1955 {
1956 strcat(fileName, ".SP3");
1957 }
1958 else if (i==3)
1959 {
1960 strcat(fileName, ".SP4");
1961 }
1962 else if (i==4)
1963 {
1964 strcat(fileName, ".SP5");
1965 }
1966 else if (i==5)
1967 {
1968 strcat(fileName, ".SP6");
1969 }
1970 else if (i==6)
1971 {
1972 strcat(fileName, ".SP7");
1973 }
1974 else if (i==7)
1975 {
1976 strcat(fileName, ".SP8");
1977 }
1978 else if (i==8)
1979 {
1980 strcat(fileName, ".SP9");
1981 }
1982 else if (i==9)
1983 {
1984 strcat(fileName, ".SPA");
1985 }
1986 else
1987 {
1988 strcat(fileName, ".SPB");
1989 }
1990 #ifdef _WIN32
1991 ifstream in(fileName , ios::binary);
1992 #else
1993 ifstream in(fileName);
1994 #endif
1995
1996 int numberOfVariables=0;
1997 if (in) // file exists
1998 {
1999 in.clear();
2000 in.seekg( 1024, ios::beg );
2001 in.read( (char*)&nextRecord,sizeof(int) );
2002 this->SwapInt(nextRecord);
2003 in.read( (char*)&numberOfRecords,sizeof(int) );
2004 this->SwapInt(numberOfRecords);
2005
2006 switch (i+1)
2007 {
2008 case 1:
2009 {
2010 numberOfVariables = 1;
2011 break;
2012 }
2013 case 2:
2014 {
2015 numberOfVariables = 2;
2016 break;
2017 }
2018 case 3:
2019 {
2020 numberOfVariables = 4;
2021 break;
2022 }
2023 case 4:
2024 {
2025 numberOfVariables = 4*this->MMAX;
2026 break;
2027 }
2028 case 5:
2029 {
2030 numberOfVariables = this->MMAX;
2031 break;
2032 }
2033 case 6:
2034 {
2035 if (this->VersionNumber <= 1.15)
2036 {
2037 numberOfVariables = 3;
2038 }
2039 else
2040 {
2041 numberOfVariables = this->MMAX + 1;
2042 }
2043 break;
2044 }
2045 case 7:
2046 {
2047 numberOfVariables = this->NMax->GetValue(0);
2048 for (int m=1; m<=this->MMAX; ++m)
2049 {
2050 numberOfVariables += this->NMax->GetValue(m);
2051 }
2052 break;
2053 }
2054 case 8:
2055 {
2056 numberOfVariables = this->MMAX;
2057 break;
2058 }
2059 case 9:
2060 {
2061 numberOfVariables = this->NumberOfScalars;
2062 break;
2063 }
2064 case 10:
2065 {
2066 numberOfVariables = this->NumberOfReactionRates;
2067 break;
2068 }
2069 case 11:
2070 {
2071 if (this->BkEpsilon)
2072 {
2073 numberOfVariables = 2;
2074 }
2075 break;
2076 }
2077 }
2078
2079 for(int j=0; j<numberOfVariables; j++)
2080 {
2081 this->VariableTimesteps->InsertValue(cnt,
2082 (nextRecord-4)/numberOfRecords);
2083 cnt++;
2084 }
2085 }
2086 in.close();
2087 }
2088 }
2089
2090 //----------------------------------------------------------------------------
MakeTimeStepTable(int numberOfVariables)2091 void vtkMFIXReader::MakeTimeStepTable(int numberOfVariables)
2092 {
2093 this->VariableTimestepTable->SetNumberOfComponents(numberOfVariables);
2094
2095 for(int i=0; i<numberOfVariables; i++)
2096 {
2097 int timestepIncrement = (int)((float)this->MaximumTimestep/
2098 (float)this->VariableTimesteps->GetValue(i) + 0.5);
2099 int timestep = 1;
2100 for (int j=0; j<this->MaximumTimestep; j++)
2101 {
2102 this->VariableTimestepTable->InsertComponent(j, i, timestep);
2103 timestepIncrement--;
2104 if (timestepIncrement <= 0)
2105 {
2106 timestepIncrement = (int)((float)this->MaximumTimestep/
2107 (float)this->VariableTimesteps->GetValue(i) + 0.5);
2108 timestep++;
2109 }
2110 if (timestep > this->VariableTimesteps->GetValue(i))
2111 {
2112 timestep = this->VariableTimesteps->GetValue(i);
2113 }
2114 }
2115 }
2116 }
2117
2118 //----------------------------------------------------------------------------
GetVariableAtTimestep(int vari,int tstep,vtkFloatArray * v)2119 void vtkMFIXReader::GetVariableAtTimestep(int vari , int tstep,
2120 vtkFloatArray *v)
2121 {
2122 // This routine opens and closes the file for each request.
2123 // Maybe keep all SPX files open, and just perform relative
2124 // moves to get to the correct location in the file
2125 // get filename that vaiable # vari is located in
2126 // assumptions : there are <10 solid phases,
2127 // <10 scalars and <10 ReactionRates (need to change this)
2128
2129 char variableName[256];
2130 strcpy(variableName, this->VariableNames->GetValue(vari));
2131 int spx = this->VariableIndexToSPX->GetValue(vari);
2132 char fileName[256];
2133
2134 for(int k=0;k<(int)sizeof(fileName);k++)
2135 {
2136 fileName[k]=0;
2137 }
2138
2139 strncpy(fileName, this->FileName, strlen(this->FileName)-4);
2140
2141 if (spx==1)
2142 {
2143 strcat(fileName, ".SP1");
2144 }
2145 else if (spx==2)
2146 {
2147 strcat(fileName, ".SP2");
2148 }
2149 else if (spx==3)
2150 {
2151 strcat(fileName, ".SP3");
2152 }
2153 else if (spx==4)
2154 {
2155 strcat(fileName, ".SP4");
2156 }
2157 else if (spx==5)
2158 {
2159 strcat(fileName, ".SP5");
2160 }
2161 else if (spx==6)
2162 {
2163 strcat(fileName, ".SP6");
2164 }
2165 else if (spx==7)
2166 {
2167 strcat(fileName, ".SP7");
2168 }
2169 else if (spx==8)
2170 {
2171 strcat(fileName, ".SP8");
2172 }
2173 else if (spx==9)
2174 {
2175 strcat(fileName, ".SP9");
2176 }
2177 else if (spx==10)
2178 {
2179 strcat(fileName, ".SPA");
2180 }
2181 else
2182 {
2183 strcat(fileName, ".SPB");
2184 }
2185
2186 int index = (vari*this->MaximumTimestep) + tstep;
2187 int nBytesSkip = this->SPXTimestepIndexTable->GetValue(index);
2188 #ifdef _WIN32
2189 ifstream in(fileName,ios::binary);
2190 #else
2191 ifstream in(fileName);
2192 #endif
2193 in.seekg(nBytesSkip,ios::beg);
2194 this->GetBlockOfFloats (in, v, this->IJKMaximum2);
2195 in.close();
2196 }
2197
2198 //----------------------------------------------------------------------------
MakeSPXTimeStepIndexTable(int nvars)2199 void vtkMFIXReader::MakeSPXTimeStepIndexTable(int nvars)
2200 {
2201
2202 int timestep;
2203 int spx;
2204 int NumberOfVariablesInSPX;
2205
2206 for (int i=0; i<nvars; i++)
2207 {
2208 for (int j=0; j<this->MaximumTimestep; j++)
2209 {
2210 timestep = (int) this->VariableTimestepTable->GetComponent(j, i);
2211 spx = this->VariableIndexToSPX->GetValue(i);
2212 NumberOfVariablesInSPX = this->SPXToNVarTable->GetValue(spx);
2213 int skip = this->VariableToSkipTable->GetValue(i);
2214 int index = (3*512) + (timestep-1) *
2215 ((NumberOfVariablesInSPX*this->SPXRecordsPerTimestep*512)+512) +
2216 512 + (skip*this->SPXRecordsPerTimestep*512);
2217 int ind = (i*this->MaximumTimestep) + j;
2218 SPXTimestepIndexTable->InsertValue(ind, index);
2219 }
2220 }
2221 }
2222
2223 //----------------------------------------------------------------------------
CalculateMaxTimeStep()2224 void vtkMFIXReader::CalculateMaxTimeStep()
2225 {
2226 this->MaximumTimestep = 0;
2227 for ( int i=0; i <= this->VariableNames->GetMaxId(); i++ )
2228 {
2229 if (this->VariableTimesteps->GetValue(i) > this->MaximumTimestep)
2230 {
2231 this->MaximumTimestep = this->VariableTimesteps->GetValue(i);
2232 }
2233 }
2234 }
2235
2236 //----------------------------------------------------------------------------
GetNumberOfVariablesInSPXFiles()2237 void vtkMFIXReader::GetNumberOfVariablesInSPXFiles()
2238 {
2239 int NumberOfVariablesInSPX = 0;
2240 int skip = 0;
2241 for (int j=1; j<this->NumberOfSPXFilesUsed; j++)
2242 {
2243 for(int i=0;i<=this->VariableNames->GetMaxId();i++)
2244 {
2245 if ((this->VariableIndexToSPX->GetValue(i) == j)
2246 && (this->VariableComponents->GetValue(i) == 1))
2247 {
2248 NumberOfVariablesInSPX++;
2249 this->VariableToSkipTable->InsertValue(i,skip);
2250 skip++;
2251 }
2252 }
2253 this->SPXToNVarTable->InsertValue(j, NumberOfVariablesInSPX);
2254 NumberOfVariablesInSPX = 0;
2255 skip = 0;
2256 }
2257 }
2258
2259 //----------------------------------------------------------------------------
FillVectorVariable(int xindex,int yindex,int zindex,vtkFloatArray * v)2260 void vtkMFIXReader::FillVectorVariable( int xindex, int yindex,
2261 int zindex, vtkFloatArray *v)
2262 {
2263 for(int i=0;i<=this->CellDataArray[xindex]->GetMaxId();i++)
2264 {
2265 v->InsertComponent(i, 0, this->CellDataArray[xindex]->GetValue(i));
2266 v->InsertComponent(i, 1, this->CellDataArray[yindex]->GetValue(i));
2267 v->InsertComponent(i, 2, this->CellDataArray[zindex]->GetValue(i));
2268 }
2269 v->Modified();
2270 }
2271
2272 //----------------------------------------------------------------------------
ConvertVectorFromCylindricalToCartesian(int xindex,int zindex)2273 void vtkMFIXReader::ConvertVectorFromCylindricalToCartesian( int xindex,
2274 int zindex)
2275 {
2276 int count = 0;
2277 double radius = 0.0;
2278 double y = 0.0;
2279 double theta = 0.0;
2280 int cnt=0;
2281
2282 for (int k=0; k< this->KMaximum2; k++)
2283 {
2284 for (int j=0; j< this->JMaximum2; j++)
2285 {
2286 for (int i=0; i< this->IMaximum2; i++)
2287 {
2288 if ( this->Flag->GetValue(cnt) < 10 )
2289 {
2290 double ucart = (this->CellDataArray[xindex]->
2291 GetValue(count)*cos(theta)) -
2292 (this->CellDataArray[zindex]->GetValue(count)*sin(theta));
2293 double wcart = (this->CellDataArray[xindex]->
2294 GetValue(count)*sin(theta)) +
2295 (this->CellDataArray[zindex]->GetValue(count)*cos(theta));
2296 this->CellDataArray[xindex]->InsertValue(count, (float)ucart);
2297 this->CellDataArray[zindex]->InsertValue(count, (float)wcart);
2298 count++;
2299 }
2300 cnt++;
2301 radius = radius + this->Dx->GetValue(i);
2302 }
2303 radius = 0.0;
2304 y = y + this->Dy->GetValue(j);
2305 }
2306 y = 0.0;
2307 theta = theta + this->Dz->GetValue(k);
2308 }
2309 }
2310
2311 //----------------------------------------------------------------------------
GetAllTimes(vtkInformationVector * outputVector)2312 void vtkMFIXReader::GetAllTimes(vtkInformationVector *outputVector)
2313 {
2314 int max = 0;
2315 int maxVar = 0;
2316
2317 for(int j=0; j<=this->VariableNames->GetMaxId(); j++)
2318 {
2319 int n = this->VariableTimesteps->GetValue(j);
2320 if (n > max)
2321 {
2322 max = n;
2323 maxVar = j;
2324 }
2325 }
2326
2327 char fileName[256];
2328
2329 for(int k=0;k<(int)sizeof(fileName);k++)
2330 {
2331 fileName[k]=0;
2332 }
2333 strncpy(fileName, this->FileName, strlen(this->FileName)-4);
2334
2335 if (maxVar==0)
2336 {
2337 strcat(fileName, ".SP1");
2338 }
2339 else if (maxVar==1)
2340 {
2341 strcat(fileName, ".SP2");
2342 }
2343 else if (maxVar==2)
2344 {
2345 strcat(fileName, ".SP3");
2346 }
2347 else if (maxVar==3)
2348 {
2349 strcat(fileName, ".SP4");
2350 }
2351 else if (maxVar==4)
2352 {
2353 strcat(fileName, ".SP5");
2354 }
2355 else if (maxVar==5)
2356 {
2357 strcat(fileName, ".SP6");
2358 }
2359 else if (maxVar==6)
2360 {
2361 strcat(fileName, ".SP7");
2362 }
2363 else if (maxVar==7)
2364 {
2365 strcat(fileName, ".SP8");
2366 }
2367 else if (maxVar==8)
2368 {
2369 strcat(fileName, ".SP9");
2370 }
2371 else if (maxVar==9)
2372 {
2373 strcat(fileName, ".SPA");
2374 }
2375 else
2376 {
2377 strcat(fileName, ".SPB");
2378 }
2379
2380 #ifdef _WIN32
2381 ifstream tfile(fileName, ios::binary);
2382 #else
2383 ifstream tfile(fileName);
2384 #endif
2385
2386 int numberOfVariablesInSPX =
2387 this->SPXToNVarTable->GetValue(this->VariableIndexToSPX->GetValue(maxVar));
2388 int offset = 512-(int)sizeof(float) +
2389 512*(numberOfVariablesInSPX*SPXRecordsPerTimestep);
2390 tfile.clear();
2391 tfile.seekg( 3*512, ios::beg ); // first time
2392 float time;
2393 double* steps = new double[this->NumberOfTimeSteps];
2394
2395 for (int i = 0; i < this->NumberOfTimeSteps; i++)
2396 {
2397 tfile.read( (char*)&time,sizeof(float) );
2398 SwapFloat(time);
2399 steps[i] = (double)time;
2400 tfile.seekg(offset,ios::cur);
2401 }
2402
2403 vtkInformation* outInfo = outputVector->GetInformationObject(0);
2404 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_STEPS(),
2405 steps, this->NumberOfTimeSteps);
2406 double timeRange[2];
2407 timeRange[0] = steps[0];
2408 timeRange[1] = steps[this->NumberOfTimeSteps - 1];
2409 outInfo->Set(vtkStreamingDemandDrivenPipeline::TIME_RANGE(), timeRange, 2);
2410
2411 tfile.close();
2412 delete [] steps;
2413 }
2414