1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    TestCONVERGECFDReader.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 // Test of vtkCONVERGECFDReader
16 
17 #include "vtkCONVERGECFDReader.h"
18 #include "vtkCellData.h"
19 #include "vtkDataArraySelection.h"
20 #include "vtkDataAssembly.h"
21 #include "vtkDebugLeaks.h"
22 #include "vtkInformation.h"
23 #include "vtkLogger.h"
24 #include "vtkPartitionedDataSetCollection.h"
25 #include "vtkPointData.h"
26 #include "vtkPolyData.h"
27 #include "vtkStreamingDemandDrivenPipeline.h"
28 #include "vtkTestUtilities.h"
29 #include "vtkUnstructuredGrid.h"
30 
31 #include <cstdlib>
32 
TestCONVERGECFDReader(int argc,char * argv[])33 int TestCONVERGECFDReader(int argc, char* argv[])
34 {
35   vtkNew<vtkCONVERGECFDReader> reader;
36 
37   {
38     // Read CONVERGE 3.1 file name
39     char* fname = vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/converge3.1-format.h5");
40     reader->SetFileName(fname);
41     reader->Update();
42     delete[] fname;
43 
44     // Check on the structure of the output partitioned dataset collection's assembly
45     auto pdc = reader->GetOutput();
46     auto assembly = pdc->GetDataAssembly();
47     if (assembly->GetNumberOfChildren(0) != 3)
48     {
49       vtkLog(ERROR, "Invalid number of streams in file.");
50       return EXIT_FAILURE;
51     }
52 
53     int stream0NodeId = assembly->GetChild(0, 0);
54     int meshNodeId = assembly->GetChild(stream0NodeId, 0);
55     int meshId = assembly->GetDataSetIndices(meshNodeId, false)[0];
56     auto mesh = vtkUnstructuredGrid::SafeDownCast(pdc->GetPartition(meshId, 0));
57     if (!mesh)
58     {
59       vtkLog(ERROR, "No mesh block found in stream 0.");
60       return EXIT_FAILURE;
61     }
62 
63     if (mesh->GetNumberOfPoints() != 176840)
64     {
65       vtkLog(ERROR,
66         "Incorrect number of points in mesh. Should be 176840, got " << mesh->GetNumberOfPoints());
67       return EXIT_FAILURE;
68     }
69 
70     if (mesh->GetNumberOfCells() != 22016)
71     {
72       vtkLog(ERROR,
73         "Incorrect number of cells in mesh. Should be 22016, got " << mesh->GetNumberOfCells());
74       return EXIT_FAILURE;
75     }
76 
77     std::vector<std::string> cellArrays = { "ASPECT_RATIO", "EPS", "EQUIV_RATIO", "FACE_WARPAGE",
78       "NON-ORTHOGONALITY", "NUM_CARTESIAN_NBRS", "NUM_INLAID_NBRS", "PRESSURE", "RANK", "SKEWNESS",
79       "STRETCH_RATIO", "TEMP_SGS", "TEMPERATURE", "TKE", "VEL_SGS", "VELOCITY" };
80 
81     if (mesh->GetCellData()->GetNumberOfArrays() != static_cast<int>(cellArrays.size()))
82     {
83       vtkLog(ERROR, "Incorrect number of cell data arrays on mesh");
84       return EXIT_FAILURE;
85     }
86 
87     for (const auto& cellArrayName : cellArrays)
88     {
89       auto cellData = mesh->GetCellData();
90       if (!cellData->HasArray(cellArrayName.c_str()))
91       {
92         vtkLog(ERROR, "Mesh is missing expected cell data array '" << cellArrayName << "'");
93         return EXIT_FAILURE;
94       }
95     }
96 
97     int surfacesNodeId = assembly->GetChild(stream0NodeId, 1);
98 
99     std::vector<std::string> pointArrays = { "FILM_FLAG", "RADIUS", "TEMP", "VELOCITY" };
100 
101     int numBlocks = assembly->GetNumberOfChildren(surfacesNodeId);
102     if (numBlocks != 42)
103     {
104       vtkLog(ERROR, "Incorrect number of surface blocks. Should be 42, got " << numBlocks);
105       return EXIT_FAILURE;
106     }
107     // Just check the first 5 surface block number of points and number of cells
108     std::string expectedBlockNames[] = { "PISTON1", "LINER1", "HEAD1", "SPARK PLUG1",
109       "SPARK PLUG ELECTRODE1" };
110     int expectedNumPoints[] = { 10095, 4159, 20202, 858, 10 };
111     int expectedNumCells[] = { 11763, 3994, 25182, 1080, 7 };
112     for (int i = 0; i < static_cast<int>(sizeof(expectedNumPoints) / sizeof(int)); ++i)
113     {
114       int surfaceNodeId = assembly->GetChild(surfacesNodeId, i);
115       int surfaceId = assembly->GetDataSetIndices(surfaceNodeId, false)[0];
116       std::string blockName(
117         pdc->GetMetaData(static_cast<unsigned int>(surfaceId))->Get(vtkCompositeDataSet::NAME()));
118       if (blockName != expectedBlockNames[i])
119       {
120         vtkLog(ERROR,
121           "Surface data block expected to be " << expectedBlockNames[i] << ", but was "
122                                                << blockName);
123         return EXIT_FAILURE;
124       }
125 
126       vtkPolyData* surface = vtkPolyData::SafeDownCast(pdc->GetPartition(surfaceId, 0));
127       if (!surface)
128       {
129         vtkLog(ERROR, "No polydata surface at block " << i);
130         return EXIT_FAILURE;
131       }
132       if (surface->GetNumberOfPoints() != expectedNumPoints[i])
133       {
134         vtkLog(ERROR, "Incorrect number of points in surface block " << i);
135         return EXIT_FAILURE;
136       }
137       if (surface->GetNumberOfCells() != expectedNumCells[i])
138       {
139         vtkLog(ERROR, "Incorrect number of cells in surface block " << i);
140         return EXIT_FAILURE;
141       }
142       if (surface->GetCellData()->GetNumberOfArrays() != static_cast<int>(cellArrays.size()))
143       {
144         vtkLog(ERROR, "Incorrect number of cell data arrays on surface at block " << i);
145         return EXIT_FAILURE;
146       }
147       for (const auto& cellArrayName : cellArrays)
148       {
149         auto cellData = surface->GetCellData();
150         if (!cellData->HasArray(cellArrayName.c_str()))
151         {
152           vtkLog(ERROR,
153             "surface " << i << " is missing expected cell data array '" << cellArrayName << "'");
154           return EXIT_FAILURE;
155         }
156       }
157     }
158 
159     int parcelNodeId = assembly->GetChild(stream0NodeId, 2);
160     int liquidParcelNodeId = assembly->GetChild(parcelNodeId, 0);
161     int liqParcel1NodeId = assembly->GetChild(liquidParcelNodeId, 0);
162     int liqParcel1Id = assembly->GetDataSetIndices(liqParcel1NodeId, false)[0];
163 
164     auto liqparcel_1 = vtkPolyData::SafeDownCast(pdc->GetPartition(liqParcel1Id, 0));
165     if (!liqparcel_1)
166     {
167       vtkLog(ERROR, "No liqparcel_1 parcel data");
168       return EXIT_FAILURE;
169     }
170 
171     if (liqparcel_1->GetNumberOfPoints() != 1581)
172     {
173       vtkLog(ERROR,
174         "Incorrect number of points in parcels. Should be 1581, got "
175           << liqparcel_1->GetNumberOfPoints());
176       return EXIT_FAILURE;
177     }
178 
179     if (liqparcel_1->GetNumberOfCells() != 1581)
180     {
181       vtkLog(ERROR,
182         "Incorrect number of cells in parcels. Should be 1581, got "
183           << liqparcel_1->GetNumberOfCells());
184       return EXIT_FAILURE;
185     }
186 
187     // Check parcel type
188     std::string expectedBlockName("LIQPARCEL_1");
189     std::string blockName(pdc->GetMetaData(liqParcel1Id)->Get(vtkCompositeDataSet::NAME()));
190     if (blockName != expectedBlockName)
191     {
192       vtkLog(
193         ERROR, "Expected block name '" << expectedBlockName << "' but got '" << blockName << "'");
194       return EXIT_FAILURE;
195     }
196 
197     // Check second stream
198     int stream1NodeId = assembly->GetChild(0, 1);
199     meshNodeId = assembly->GetChild(stream1NodeId, 0);
200     meshId = assembly->GetDataSetIndices(meshNodeId, false)[0];
201     mesh = vtkUnstructuredGrid::SafeDownCast(pdc->GetPartition(meshId, 0));
202     if (!mesh)
203     {
204       vtkLog(ERROR, "No mesh block found in stream 1.");
205       return EXIT_FAILURE;
206     }
207 
208     if (mesh->GetNumberOfPoints() != 178273)
209     {
210       vtkLog(ERROR,
211         "Incorrect number of points in mesh. Should be 178273, got " << mesh->GetNumberOfPoints());
212       return EXIT_FAILURE;
213     }
214 
215     if (mesh->GetNumberOfCells() != 22369)
216     {
217       vtkLog(ERROR,
218         "Incorrect number of cells in mesh. Should be 22369, got " << mesh->GetNumberOfCells());
219       return EXIT_FAILURE;
220     }
221 
222     surfacesNodeId = assembly->GetChild(stream1NodeId, 1);
223 
224     numBlocks = assembly->GetNumberOfChildren(surfacesNodeId);
225     if (numBlocks != 42)
226     {
227       vtkLog(ERROR, "Incorrect number of surface blocks. Should be 42, got " << numBlocks);
228       return EXIT_FAILURE;
229     }
230 
231     parcelNodeId = assembly->GetChild(stream1NodeId, 2);
232     liquidParcelNodeId = assembly->GetChild(parcelNodeId, 0);
233     liqParcel1NodeId = assembly->GetChild(liquidParcelNodeId, 0);
234     liqParcel1Id = assembly->GetDataSetIndices(liqParcel1NodeId, false)[0];
235 
236     liqparcel_1 = vtkPolyData::SafeDownCast(pdc->GetPartition(liqParcel1Id, 0));
237     if (!liqparcel_1)
238     {
239       vtkLog(ERROR, "No liqparcel_1 parcel data");
240       return EXIT_FAILURE;
241     }
242 
243     if (liqparcel_1->GetNumberOfPoints() != 1798)
244     {
245       vtkLog(ERROR,
246         "Incorrect number of points in parcels. Should be 1798, got "
247           << liqparcel_1->GetNumberOfPoints());
248       return EXIT_FAILURE;
249     }
250 
251     if (liqparcel_1->GetNumberOfCells() != 1798)
252     {
253       vtkLog(ERROR,
254         "Incorrect number of cells in parcels. Should be 1798, got "
255           << liqparcel_1->GetNumberOfCells());
256       return EXIT_FAILURE;
257     }
258 
259     // Check third stream
260     int stream2NodeId = assembly->GetChild(0, 2);
261     meshNodeId = assembly->GetChild(stream2NodeId, 0);
262     meshId = assembly->GetDataSetIndices(meshNodeId, false)[0];
263     mesh = vtkUnstructuredGrid::SafeDownCast(pdc->GetPartition(meshId, 0));
264     if (!mesh)
265     {
266       vtkLog(ERROR, "No mesh block found in stream 2.");
267       return EXIT_FAILURE;
268     }
269 
270     if (mesh->GetNumberOfPoints() != 3620)
271     {
272       vtkLog(ERROR,
273         "Incorrect number of points in mesh. Should be 3620, got " << mesh->GetNumberOfPoints());
274       return EXIT_FAILURE;
275     }
276 
277     if (mesh->GetNumberOfCells() != 124)
278     {
279       vtkLog(ERROR,
280         "Incorrect number of cells in mesh. Should be 124, got " << mesh->GetNumberOfCells());
281       return EXIT_FAILURE;
282     }
283 
284     // Ensure there are no parcels in the third stream
285     if (assembly->GetNumberOfChildren(stream2NodeId) != 2)
286     {
287       vtkLog(ERROR,
288         "Number of children should be 2, but is " << assembly->GetNumberOfChildren(stream2NodeId)
289                                                   << " instead");
290       return EXIT_FAILURE;
291     }
292   }
293 
294   {
295     // Read file name (CONVERGE 3.0 file).
296     char* fname = vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/post_5016_spray.h5");
297     reader->GetCellDataArraySelection()->RemoveAllArrays();
298     reader->GetParcelDataArraySelection()->RemoveAllArrays();
299     reader->SetFileName(fname);
300     reader->Update();
301     delete[] fname;
302 
303     // Check on the structure of the output multiblock dataset
304     auto pdc = reader->GetOutput();
305     auto assembly = pdc->GetDataAssembly();
306     if (assembly->GetNumberOfChildren(0) != 1)
307     {
308       vtkLog(ERROR, "Invalid number of streams in file.");
309       return EXIT_FAILURE;
310     }
311 
312     int stream0NodeId = assembly->GetChild(0, 0);
313     int meshNodeId = assembly->GetChild(stream0NodeId, 0);
314     int meshId = assembly->GetDataSetIndices(meshNodeId, false)[0];
315     auto mesh = vtkUnstructuredGrid::SafeDownCast(pdc->GetPartition(meshId, 0));
316     if (!mesh)
317     {
318       vtkLog(ERROR, "No mesh block found in file.");
319       return EXIT_FAILURE;
320     }
321 
322     if (mesh->GetNumberOfPoints() != 12242)
323     {
324       vtkLog(ERROR, "Incorrect number of points in mesh.");
325       return EXIT_FAILURE;
326     }
327 
328     if (mesh->GetNumberOfCells() != 3378)
329     {
330       vtkLog(ERROR, "Incorrect number of cells in mesh.");
331       return EXIT_FAILURE;
332     }
333 
334     std::vector<std::string> cellArrays = { "DENSITY", "EPS", "EQUIV_RATIO", "LAMBDA",
335       "MASSFRAC_C7H16", "MASSFRAC_CO", "MASSFRAC_CO2", "MASSFRAC_H2", "MASSFRAC_H2O", "MASSFRAC_O2",
336       "PRESSURE", "RANK", "REACT_RATIO", "SIE", "TEMPERATURE", "TKE", "VELOCITY", "VISC" };
337 
338     std::vector<std::string> pointArrays = { "FILM_FLAG", "RADIUS", "TEMP", "VELOCITY" };
339 
340     if (mesh->GetCellData()->GetNumberOfArrays() != static_cast<int>(cellArrays.size()))
341     {
342       vtkLog(ERROR, "Incorrect number of cell data arrays on mesh");
343       return EXIT_FAILURE;
344     }
345 
346     for (const auto& cellArrayName : cellArrays)
347     {
348       auto cellData = mesh->GetCellData();
349       if (!cellData->HasArray(cellArrayName.c_str()))
350       {
351         vtkLog(ERROR, "Mesh is missing expected cell data array '" << cellArrayName << "'");
352         return EXIT_FAILURE;
353       }
354     }
355 
356     int surfacesNodeId = assembly->GetChild(stream0NodeId, 1);
357     int numBlocks = assembly->GetNumberOfChildren(surfacesNodeId);
358     if (numBlocks != 7)
359     {
360       vtkLog(ERROR, "Incorrect number of surface blocks. Should be 7, got " << numBlocks);
361       return EXIT_FAILURE;
362     }
363     int expectedNumPoints[] = { 5535, 837, 829, 510, 1374, 0, 0 };
364     int expectedNumCells[] = { 6038, 770, 763, 461, 1286, 0, 0 };
365     for (int i = 0; i < numBlocks; ++i)
366     {
367       int surfaceNodeId = assembly->GetChild(surfacesNodeId, i);
368       int surfaceId = assembly->GetDataSetIndices(surfaceNodeId, false)[0];
369       vtkPolyData* surface = vtkPolyData::SafeDownCast(pdc->GetPartition(surfaceId, 0));
370       if (!surface)
371       {
372         vtkLog(ERROR, "No polydata surface at block " << i);
373         return EXIT_FAILURE;
374       }
375       if (surface->GetNumberOfPoints() != expectedNumPoints[i])
376       {
377         vtkLog(ERROR, "Incorrect number of points in surface block " << i);
378         return EXIT_FAILURE;
379       }
380       if (surface->GetNumberOfCells() != expectedNumCells[i])
381       {
382         vtkLog(ERROR, "Incorrect number of cells in surface block " << i);
383         return EXIT_FAILURE;
384       }
385       if (surface->GetCellData()->GetNumberOfArrays() != static_cast<int>(cellArrays.size()))
386       {
387         vtkLog(ERROR, "Incorrect number of cell data arrays on surface at block " << i);
388         return EXIT_FAILURE;
389       }
390       for (const auto& cellArrayName : cellArrays)
391       {
392         auto cellData = surface->GetCellData();
393         if (!cellData->HasArray(cellArrayName.c_str()))
394         {
395           vtkLog(ERROR,
396             "surface " << i << " is missing expected cell data array '" << cellArrayName << "'");
397           return EXIT_FAILURE;
398         }
399       }
400     }
401 
402     int parcelsNodeId = assembly->GetChild(stream0NodeId, 2);
403     int parcelsId = assembly->GetDataSetIndices(parcelsNodeId, false)[0];
404     auto parcels = vtkPolyData::SafeDownCast(pdc->GetPartition(parcelsId, 0));
405     if (!parcels)
406     {
407       vtkLog(ERROR, "No parcels block found in file.");
408       return EXIT_FAILURE;
409     }
410 
411     if (parcels->GetNumberOfPoints() != 185732)
412     {
413       vtkLog(ERROR, "Incorrect number of points in parcels.");
414       return EXIT_FAILURE;
415     }
416 
417     if (parcels->GetNumberOfCells() != 185732)
418     {
419       vtkLog(ERROR, "Incorrect number of cells in parcels.");
420       return EXIT_FAILURE;
421     }
422 
423     if (parcels->GetPointData()->GetNumberOfArrays() != static_cast<int>(pointArrays.size()))
424     {
425       vtkLog(ERROR, "Incorrect number of parcel data arrays");
426       return EXIT_FAILURE;
427     }
428     for (const auto& pointArrayName : pointArrays)
429     {
430       auto pointData = parcels->GetPointData();
431       if (!pointData->HasArray(pointArrayName.c_str()))
432       {
433         vtkLog(ERROR, "Parcels are missing expected point data array '" << pointArrayName << "'");
434         return EXIT_FAILURE;
435       }
436     }
437 
438     // Test array selection
439     auto cellArraySelection = reader->GetCellDataArraySelection();
440     cellArraySelection->DisableArray("EPS");
441     cellArraySelection->DisableArray("DENSITY");
442     auto parcelArraySelection = reader->GetParcelDataArraySelection();
443     parcelArraySelection->DisableArray("RADIUS");
444     reader->Update();
445 
446     pdc = reader->GetOutput();
447     assembly = pdc->GetDataAssembly();
448     mesh = vtkUnstructuredGrid::SafeDownCast(pdc->GetPartition(meshId, 0));
449     parcels = vtkPolyData::SafeDownCast(pdc->GetPartition(parcelsId, 0));
450 
451     auto cellData = mesh->GetCellData();
452     if (cellData->GetArray("EPS") != nullptr)
453     {
454       vtkLog(ERROR, "Data array 'EPS' should not have been read but is available");
455       return EXIT_FAILURE;
456     }
457     if (cellData->GetArray("DENSITY") != nullptr)
458     {
459       vtkLog(ERROR, "Data array 'DENSITY' should not have been read but is available");
460       return EXIT_FAILURE;
461     }
462 
463     for (int i = 0; i < numBlocks; ++i)
464     {
465       int surfaceNodeId = assembly->GetChild(surfacesNodeId, i);
466       int surfaceId = assembly->GetDataSetIndices(surfaceNodeId, false)[0];
467       vtkPolyData* surface = vtkPolyData::SafeDownCast(pdc->GetPartition(surfaceId, 0));
468       if (!surface)
469       {
470         vtkLog(ERROR, "No polydata surface at block " << i);
471         return EXIT_FAILURE;
472       }
473 
474       cellData = surface->GetCellData();
475       if (cellData->GetArray("EPS") != nullptr)
476       {
477         vtkLog(ERROR, "Data array 'EPS' should not have been read but is available");
478         return EXIT_FAILURE;
479       }
480       if (cellData->GetArray("DENSITY") != nullptr)
481       {
482         vtkLog(ERROR, "Data array 'DENSITY' should not have been read but is available");
483         return EXIT_FAILURE;
484       }
485     }
486 
487     auto parcelData = parcels->GetPointData();
488     if (parcelData->GetArray("RADIUS") != nullptr)
489     {
490       vtkLog(ERROR, "Data array 'RADIUS' should not have been read but is available");
491       return EXIT_FAILURE;
492     }
493   }
494 
495   return EXIT_SUCCESS;
496 }
497