1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    TestTRUCHASReader.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 // Input test to validate ability to read GE TRUCHAS files
17 #include "vtkCellData.h"
18 #include "vtkDoubleArray.h"
19 #include "vtkMultiBlockDataSet.h"
20 #include "vtkNew.h"
21 #include "vtkPointData.h"
22 #include "vtkTRUCHASReader.h"
23 #include "vtkTestUtilities.h"
24 #include "vtkTesting.h"
25 #include "vtkUnstructuredGrid.h"
26 
27 #include "vtkInformation.h"
28 #include "vtkStreamingDemandDrivenPipeline.h"
29 
30 namespace
31 {
AE(double v1,double v2)32 bool AE(double v1, double v2)
33 {
34   if (fabs(v2 - v1) > 0.001)
35   {
36     cerr << v2 << "!=" << v1 << endl;
37     return false;
38   }
39   return true;
40 }
41 }
42 
TestTRUCHASReader(int argc,char * argv[])43 int TestTRUCHASReader(int argc, char* argv[])
44 {
45   vtkNew<vtkTesting> testing;
46   testing->AddArguments(argc, argv);
47   char* fileName =
48     vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/TRUCHAS/viscoplastic-ring.h5");
49 
50   vtkNew<vtkTRUCHASReader> reader;
51   reader->SetFileName(fileName);
52   delete[] fileName;
53   reader->UpdateInformation();
54   int nb = reader->GetNumberOfBlockArrays();
55   cerr << nb << " BLOCKS" << endl;
56   for (int b = 0; b < nb; b++)
57   {
58     cerr << "BLOCK ID " << b << " named " << reader->GetBlockArrayName(b) << endl;
59   }
60   reader->SetBlockArrayStatus("1", 0); // block nums start at 1
61   reader->SetBlockArrayStatus("2", 1); // block nums start at 1
62   reader->SetBlockArrayStatus("3", 0); // block nums start at 1
63 
64   int nca = reader->GetNumberOfCellArrays();
65   cerr << nca << " CELL ARRAYS" << endl;
66   for (int a = 0; a < nca; a++)
67   {
68     cerr << "ARRAY " << a << " named " << reader->GetCellArrayName(a) << endl;
69   }
70   cerr << "IGNORE VOF" << endl;
71   reader->SetCellArrayStatus("VOF", 0);
72 
73   int npa = reader->GetNumberOfPointArrays();
74   cerr << npa << " POINT ARRAYS" << endl;
75   for (int a = 0; a < npa; a++)
76   {
77     cerr << "ARRAY " << a << " named " << reader->GetPointArrayName(a) << endl;
78   }
79   cerr << "IGNORE Displacement" << endl;
80   reader->SetPointArrayStatus("Displacement", 0);
81   reader->Update();
82 
83   vtkUnstructuredGrid* grid = vtkUnstructuredGrid::SafeDownCast(reader->GetOutput()->GetBlock(1));
84   if (!grid)
85   {
86     cerr << "Could not open first block of known good file" << endl;
87     return EXIT_FAILURE;
88   }
89   int rnb = 0; // we produce empty blocks when deselected
90   for (int b = 0; b < nb; b++)
91   {
92     if (reader->GetOutput()->GetBlock(b))
93     {
94       rnb++;
95     }
96   }
97   if ((rnb != 1) || (reader->GetOutput()->GetNumberOfBlocks() != 3))
98   {
99     cerr << "Got unexpected number of blocks, found " << rnb << "/"
100          << reader->GetOutput()->GetNumberOfBlocks() << " instead of " << 1 << "/" << 3 << endl;
101   }
102 
103   cerr << "---- CELL ARRAYS ----" << endl;
104   const int expectedNumCArrays = nca - 1;
105   if (nca > 0 && grid->GetCellData()->GetNumberOfArrays() != expectedNumCArrays)
106   {
107     cerr << "Got unexpected number of cell arrays, found "
108          << grid->GetCellData()->GetNumberOfArrays() << " instead of " << expectedNumCArrays
109          << endl;
110     return EXIT_FAILURE;
111   }
112   for (int a = 0; a < grid->GetCellData()->GetNumberOfArrays(); a++)
113   {
114     vtkDataArray* da = grid->GetCellData()->GetArray(a);
115     cerr << da->GetName() << endl;
116   }
117 
118   cerr << "---- POINT ARRAYS ----" << endl;
119   const int expectedNumPArrays = npa - 1;
120   if (npa > 0 && grid->GetPointData()->GetNumberOfArrays() != expectedNumPArrays)
121   {
122     cerr << "Got unexpected number of point arrays, found "
123          << grid->GetPointData()->GetNumberOfArrays() << " instead of " << expectedNumPArrays
124          << endl;
125     return EXIT_FAILURE;
126   }
127   for (int a = 0; a < grid->GetPointData()->GetNumberOfArrays(); a++)
128   {
129     vtkDataArray* da = grid->GetPointData()->GetArray(a);
130     cerr << da->GetName() << endl;
131   }
132 
133   const int expectedNumPoints = 496;
134   if (grid->GetNumberOfPoints() != expectedNumPoints)
135   {
136     cerr << "Got unexpected number of points from file " << grid->GetNumberOfPoints()
137          << " instead of " << expectedNumPoints << endl;
138     return EXIT_FAILURE;
139   }
140   const int expectedNumCells = 180;
141   if (grid->GetNumberOfCells() != expectedNumCells)
142   {
143     cerr << "Got unexpected number of cells from file " << grid->GetNumberOfCells()
144          << " instead of " << expectedNumCells << endl;
145     return EXIT_FAILURE;
146   }
147 
148   vtkDoubleArray* da = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("Grad_T"));
149   if (!da)
150   {
151     cerr << "Couldn't get "
152          << "Grad_T"
153          << " array" << endl;
154     return EXIT_FAILURE;
155   }
156 
157   double* ptr = da->GetTuple(42);
158   const double eVals[3] = { -10.4436, -4.32586, -10.4913 };
159 
160   if (!AE(*ptr, eVals[0]) || !AE(*(ptr + 1), eVals[1]) || !AE(*(ptr + 2), eVals[2]))
161   {
162     cerr << "Got unexpected values from Grad_T array for cell 42 " << *ptr << "," << *(ptr + 1)
163          << "," << *(ptr + 2) << " instead of " << eVals[0] << "," << eVals[1] << "," << eVals[2]
164          << endl;
165     return EXIT_FAILURE;
166   }
167   reader->SetCellArrayStatus("VOF", 1);
168   reader->Update();
169   if (grid->GetCellData()->GetNumberOfArrays() != expectedNumCArrays + 1)
170   {
171     cerr << "Got unexpected number of cell arrays, found "
172          << grid->GetCellData()->GetNumberOfArrays() << " instead of " << expectedNumCArrays + 1
173          << endl;
174     return EXIT_FAILURE;
175   }
176 
177   vtkInformation* inf = reader->GetExecutive()->GetOutputInformation(0);
178   int numTimes = inf->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS());
179   cerr << "FOUND " << numTimes << " timesteps " << endl;
180   double tAlpha = inf->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), 0);
181   double tOmega = inf->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), numTimes - 1);
182   const int expectedNumTimes = 2;
183   const double expectedMinT = 0.0;
184   const double expectedMaxT = 0.5;
185   if (numTimes != 2 || !AE(tAlpha, expectedMinT) || !AE(tOmega, expectedMaxT))
186   {
187     cerr << "Got unexpected times." << endl;
188     cerr << numTimes << " not " << expectedNumTimes << " times ";
189     cerr << tAlpha << " not " << expectedMinT << " first time";
190     cerr << tOmega << " not " << expectedMaxT << " last time";
191     return EXIT_FAILURE;
192   }
193   const int divs = 3;
194   const double expectedRanges[divs][2] = {
195     { 0, 0 },               // before
196     { 0, 0 },               // after first
197     { -1.99025, -0.85729 }, // after second
198   };
199 
200   for (int i = 0; i < divs; i++)
201   {
202     double tNext = tAlpha - 0.1 + i * (tOmega - tAlpha) * 2.5 / divs;
203     reader->UpdateTimeStep(tNext);
204     grid = vtkUnstructuredGrid::SafeDownCast(reader->GetOutput()->GetBlock(1));
205     da = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("dTdt"));
206     double* mM = da->GetRange();
207     cerr << "ts " << i << ":" << tNext << " got " << mM[0] << "," << mM[1] << endl;
208     if (!AE(mM[0], expectedRanges[i][0]) || !AE(mM[1], expectedRanges[i][1]))
209     {
210       cerr << "Got unexpected ranges at time " << tNext << " " << mM[0] << "," << mM[1]
211            << " instead of " << expectedRanges[i][0] << "," << expectedRanges[i][1] << endl;
212       return EXIT_FAILURE;
213     }
214   }
215 
216   return EXIT_SUCCESS;
217 }
218