1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkXdmfHeavyData.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 #include "vtkXdmfHeavyData.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkCellTypes.h"
20 #include "vtkDataArrayRange.h"
21 #include "vtkDataObjectTypes.h"
22 #include "vtkDoubleArray.h"
23 #include "vtkExtractSelectedIds.h"
24 #include "vtkFloatArray.h"
25 #include "vtkIdTypeArray.h"
26 #include "vtkInformation.h"
27 #include "vtkMath.h"
28 #include "vtkMergePoints.h"
29 #include "vtkMultiBlockDataSet.h"
30 #include "vtkNew.h"
31 #include "vtkObjectFactory.h"
32 #include "vtkPointData.h"
33 #include "vtkPolyData.h"
34 #include "vtkRectilinearGrid.h"
35 #include "vtkSelection.h"
36 #include "vtkSelectionNode.h"
37 #include "vtkSmartPointer.h"
38 #include "vtkStructuredData.h"
39 #include "vtkStructuredGrid.h"
40 #include "vtkUniformGrid.h"
41 #include "vtkUnsignedCharArray.h"
42 #include "vtkUnstructuredGrid.h"
43 #include "vtkXdmfDataArray.h"
44 #include "vtkXdmfReader.h"
45 #include "vtkXdmfReaderInternal.h"
46
47 #include <algorithm>
48 #include <cassert>
49 #include <deque>
50 #include <numeric>
51 #include <type_traits>
52 #include <vector>
53
54 #include "vtk_libxml2.h"
55 #include VTKLIBXML2_HEADER(tree.h)
56
57 #ifdef VTK_USE_64BIT_IDS
58 typedef XdmfInt64 vtkXdmfIdType;
59 #else
60 typedef XdmfInt32 vtkXdmfIdType;
61 #endif
62
63 using namespace xdmf2;
64
vtkScaleExtents(int in_exts[6],int out_exts[6],int stride[3])65 static void vtkScaleExtents(int in_exts[6], int out_exts[6], int stride[3])
66 {
67 out_exts[0] = in_exts[0] / stride[0];
68 out_exts[1] = in_exts[1] / stride[0];
69 out_exts[2] = in_exts[2] / stride[1];
70 out_exts[3] = in_exts[3] / stride[1];
71 out_exts[4] = in_exts[4] / stride[2];
72 out_exts[5] = in_exts[5] / stride[2];
73 }
74
vtkGetDims(int exts[6],int dims[3])75 static void vtkGetDims(int exts[6], int dims[3])
76 {
77 dims[0] = exts[1] - exts[0] + 1;
78 dims[1] = exts[3] - exts[2] + 1;
79 dims[2] = exts[5] - exts[4] + 1;
80 }
81
82 //------------------------------------------------------------------------------
vtkXdmfHeavyData(vtkXdmfDomain * domain,vtkAlgorithm * reader)83 vtkXdmfHeavyData::vtkXdmfHeavyData(vtkXdmfDomain* domain, vtkAlgorithm* reader)
84 {
85 this->Reader = reader;
86 this->Piece = 0;
87 this->NumberOfPieces = 0;
88 this->GhostLevels = 0;
89 this->Extents[0] = this->Extents[2] = this->Extents[4] = 0;
90 this->Extents[1] = this->Extents[3] = this->Extents[5] = -1;
91 this->Domain = domain;
92 this->Stride[0] = this->Stride[1] = this->Stride[2] = 1;
93 }
94
95 //------------------------------------------------------------------------------
96 vtkXdmfHeavyData::~vtkXdmfHeavyData() = default;
97
98 //------------------------------------------------------------------------------
ReadData()99 vtkDataObject* vtkXdmfHeavyData::ReadData()
100 {
101 if (this->Domain->GetNumberOfGrids() == 1)
102 {
103 // There's just 1 grid. Now in serial, this is all good. In parallel, we
104 // need to be care:
105 // 1. If the data is structured, we respect the update-extent and read
106 // accordingly.
107 // 2. If the data is unstructrued, we read only on the root node. The user
108 // can apply D3 or something to repartition the data.
109 return this->ReadData(this->Domain->GetGrid(0));
110 }
111
112 // this code is similar to ReadComposite() however we cannot use the same code
113 // since the API for getting the children differs on the domain and the grid.
114
115 bool distribute_leaf_nodes = this->NumberOfPieces > 1;
116 XdmfInt32 numChildren = this->Domain->GetNumberOfGrids();
117 int number_of_leaf_nodes = 0;
118
119 vtkMultiBlockDataSet* mb = vtkMultiBlockDataSet::New();
120 mb->SetNumberOfBlocks(numChildren);
121
122 for (XdmfInt32 cc = 0; cc < numChildren; cc++)
123 {
124 XdmfGrid* xmfChild = this->Domain->GetGrid(cc);
125 mb->GetMetaData(cc)->Set(vtkCompositeDataSet::NAME(), xmfChild->GetName());
126 bool child_is_leaf = (xmfChild->IsUniform() != 0);
127 if (!child_is_leaf || !distribute_leaf_nodes ||
128 (number_of_leaf_nodes % this->NumberOfPieces) == this->Piece)
129 {
130 // it's possible that the data has way too many blocks, in which case the
131 // reader didn't present the user with capabilities to select the actual
132 // leaf node blocks as is the norm, instead only top-level grids were
133 // shown. In that case we need to ensure that we skip grids the user
134 // wanted us to skip explicitly.
135 if (!this->Domain->GetGridSelection()->ArrayIsEnabled(xmfChild->GetName()))
136 {
137 continue;
138 }
139 vtkDataObject* childDO = this->ReadData(xmfChild);
140 if (childDO)
141 {
142 mb->SetBlock(cc, childDO);
143 childDO->Delete();
144 }
145 }
146 number_of_leaf_nodes += child_is_leaf ? 1 : 0;
147 }
148
149 return mb;
150 }
151
152 //------------------------------------------------------------------------------
ReadData(XdmfGrid * xmfGrid,int blockId)153 vtkDataObject* vtkXdmfHeavyData::ReadData(XdmfGrid* xmfGrid, int blockId)
154 {
155 if (!xmfGrid || xmfGrid->GetGridType() == XDMF_GRID_UNSET)
156 {
157 // sanity check-ensure that the xmfGrid is valid.
158 return nullptr;
159 }
160
161 XdmfInt32 gridType = (xmfGrid->GetGridType() & XDMF_GRID_MASK);
162 if (gridType == XDMF_GRID_COLLECTION &&
163 xmfGrid->GetCollectionType() == XDMF_GRID_COLLECTION_TEMPORAL)
164 {
165 // grid is a temporal collection, pick the sub-grid with matching time and
166 // process that.
167 return this->ReadTemporalCollection(xmfGrid, blockId);
168 }
169 else if (gridType == XDMF_GRID_COLLECTION || gridType == XDMF_GRID_TREE)
170 {
171 return this->ReadComposite(xmfGrid);
172 }
173
174 // grid is a primitive grid, so read the data.
175 return this->ReadUniformData(xmfGrid, blockId);
176 }
177
178 //------------------------------------------------------------------------------
ReadComposite(XdmfGrid * xmfComposite)179 vtkDataObject* vtkXdmfHeavyData::ReadComposite(XdmfGrid* xmfComposite)
180 {
181 assert(((xmfComposite->GetGridType() & XDMF_GRID_COLLECTION &&
182 xmfComposite->GetCollectionType() != XDMF_GRID_COLLECTION_TEMPORAL) ||
183 (xmfComposite->GetGridType() & XDMF_GRID_TREE)) &&
184 "Input must be a spatial collection or a tree");
185
186 vtkMultiBlockDataSet* multiBlock = vtkMultiBlockDataSet::New();
187 XdmfInt32 numChildren = xmfComposite->GetNumberOfChildren();
188 multiBlock->SetNumberOfBlocks(numChildren);
189
190 bool distribute_leaf_nodes =
191 (xmfComposite->GetGridType() & XDMF_GRID_COLLECTION && this->NumberOfPieces > 1);
192
193 int number_of_leaf_nodes = 0;
194 for (XdmfInt32 cc = 0; cc < numChildren; cc++)
195 {
196 XdmfGrid* xmfChild = xmfComposite->GetChild(cc);
197 multiBlock->GetMetaData(cc)->Set(vtkCompositeDataSet::NAME(), xmfChild->GetName());
198 bool child_is_leaf = (xmfChild->IsUniform() != 0);
199 if (!child_is_leaf || !distribute_leaf_nodes ||
200 (number_of_leaf_nodes % this->NumberOfPieces) == this->Piece)
201 {
202 vtkDataObject* childDO = this->ReadData(xmfChild, cc);
203 if (childDO)
204 {
205 multiBlock->SetBlock(cc, childDO);
206 childDO->Delete();
207 }
208 }
209 number_of_leaf_nodes += child_is_leaf ? 1 : 0;
210 }
211
212 return multiBlock;
213 }
214
215 //------------------------------------------------------------------------------
ReadTemporalCollection(XdmfGrid * xmfTemporalCollection,int blockId)216 vtkDataObject* vtkXdmfHeavyData::ReadTemporalCollection(
217 XdmfGrid* xmfTemporalCollection, int blockId)
218 {
219 assert(xmfTemporalCollection->GetGridType() & XDMF_GRID_COLLECTION &&
220 xmfTemporalCollection->GetCollectionType() == XDMF_GRID_COLLECTION_TEMPORAL &&
221 "Input must be a temporal collection");
222
223 // Find the children that are valid for the requested time (this->Time) and
224 // read only those.
225
226 // FIXME: I am tempted to remove support for supporting multiple matching
227 // sub-grids for a time-step since that changes the composite data hierarchy
228 // over time which makes it hard to use filters such as vtkExtractBlock etc.
229
230 std::deque<XdmfGrid*> valid_children;
231 for (XdmfInt32 cc = 0; cc < xmfTemporalCollection->GetNumberOfChildren(); cc++)
232 {
233 XdmfGrid* child = xmfTemporalCollection->GetChild(cc);
234 if (child)
235 {
236 // ensure that we set correct epsilon for comparison
237 // BUG #0013766.
238 child->GetTime()->SetEpsilon(VTK_DBL_EPSILON);
239 if (child->GetTime()->IsValid(this->Time, this->Time))
240 {
241 valid_children.push_back(child);
242 }
243 }
244 }
245 // if no child matched this timestep, handle the case where the user didn't
246 // specify any <Time /> element for the temporal collection.
247 for (XdmfInt32 cc = 0;
248 valid_children.empty() && cc < xmfTemporalCollection->GetNumberOfChildren(); cc++)
249 {
250 XdmfGrid* child = xmfTemporalCollection->GetChild(cc);
251 if (child && child->GetTime()->GetTimeType() == XDMF_TIME_UNSET)
252 {
253 valid_children.push_back(child);
254 }
255 }
256
257 if (valid_children.empty())
258 {
259 return nullptr;
260 }
261
262 std::deque<vtkSmartPointer<vtkDataObject>> child_data_objects;
263 std::deque<XdmfGrid*>::iterator iter;
264 for (iter = valid_children.begin(); iter != valid_children.end(); ++iter)
265 {
266 vtkDataObject* childDO = this->ReadData(*iter, blockId);
267 if (childDO)
268 {
269 child_data_objects.emplace_back(childDO);
270 childDO->Delete();
271 }
272 }
273
274 if (child_data_objects.size() == 1)
275 {
276 vtkDataObject* dataObject = child_data_objects[0];
277 dataObject->Register(nullptr);
278 return dataObject;
279 }
280 else if (child_data_objects.size() > 1)
281 {
282 vtkMultiBlockDataSet* mb = vtkMultiBlockDataSet::New();
283 mb->SetNumberOfBlocks(static_cast<unsigned int>(child_data_objects.size()));
284 for (unsigned int cc = 0; cc < static_cast<unsigned int>(child_data_objects.size()); cc++)
285 {
286 mb->SetBlock(cc, child_data_objects[cc]);
287 }
288 return mb;
289 }
290
291 return nullptr;
292 }
293
294 //------------------------------------------------------------------------------
295 // Read a non-composite grid. Note here uniform has nothing to do with
296 // vtkUniformGrid but to what Xdmf's GridType="Uniform".
ReadUniformData(XdmfGrid * xmfGrid,int blockId)297 vtkDataObject* vtkXdmfHeavyData::ReadUniformData(XdmfGrid* xmfGrid, int blockId)
298 {
299 assert(xmfGrid->IsUniform() && "Input must be a uniform xdmf grid.");
300
301 int vtk_data_type = this->Domain->GetVTKDataType(xmfGrid);
302
303 if (!this->Domain->GetGridSelection()->ArrayIsEnabled(xmfGrid->GetName()))
304 {
305 // simply create an empty data-object of the correct type and return it.
306 return vtkDataObjectTypes::NewDataObject(vtk_data_type);
307 }
308
309 // Read heavy data for grid geometry/topology. This does not read any
310 // data-arrays. They are read explicitly.
311 XdmfTopology* topo = xmfGrid->GetTopology();
312 XdmfGeometry* geom = xmfGrid->GetGeometry();
313 xmlChar* filePtr;
314
315 bool caching = true;
316 XdmfDOM* topoDom = topo->GetDOM();
317 XdmfXmlNode topoNode = topo->GetElement();
318 XdmfXmlNode topoNodeDataItem = topoDom->FindElement("DataItem", 0, topoNode);
319 std::string topoFilename = "NULL";
320 if (topoNodeDataItem && caching)
321 {
322 filePtr = topoNodeDataItem->children->content;
323 if (filePtr != nullptr)
324 {
325 topoFilename = reinterpret_cast<char*>(filePtr);
326 }
327 else
328 {
329 // vtkErrorWithObjectMacro(this->Reader, << "Cannot find DataItem element in topology xml, no
330 // caching possible");
331 caching = false;
332 }
333 }
334 else
335 {
336 caching = false;
337 }
338
339 XdmfDOM* geomDom = geom->GetDOM();
340 XdmfXmlNode geomNode = geom->GetElement();
341 XdmfXmlNode geomNodeDataItem = geomDom->FindElement("DataItem", 0, geomNode);
342 std::string geomFilename = "NULL";
343 if (geomNodeDataItem && caching)
344 {
345 filePtr = geomNodeDataItem->children->content;
346 if (filePtr != nullptr)
347 {
348 geomFilename = reinterpret_cast<char*>(filePtr);
349 }
350 else
351 {
352 vtkErrorWithObjectMacro(
353 this->Reader, << "Cannot find DataItem element in geometry xml, no caching possible");
354 caching = false;
355 }
356 }
357 else
358 {
359 caching = false;
360 }
361
362 vtkXdmfReader::XdmfReaderCachedData& cache =
363 vtkXdmfReader::SafeDownCast(this->Reader)->GetDataSetCache();
364 vtkXdmfReader::XdmfDataSetTopoGeoPath& cachedData = cache[blockId];
365 if (caching && (cachedData.topologyPath == topoFilename) &&
366 (cachedData.geometryPath == geomFilename))
367 {
368 vtkDataSet* ds = vtkDataSet::SafeDownCast(
369 vtkDataObjectTypes::NewDataObject(cachedData.dataset->GetDataObjectType()));
370 ds->ShallowCopy(cachedData.dataset);
371 this->ReadAttributes(ds, xmfGrid);
372 return ds;
373 }
374
375 if (caching)
376 {
377 cachedData.topologyPath = topoFilename;
378 cachedData.geometryPath = geomFilename;
379 if (cache[blockId].dataset != nullptr)
380 {
381 cache[blockId].dataset->Delete();
382 cache[blockId].dataset = nullptr;
383 }
384 }
385
386 XdmfInt32 status = xmfGrid->Update();
387 if (status == XDMF_FAIL)
388 {
389 return nullptr;
390 }
391
392 vtkDataObject* dataObject = nullptr;
393
394 switch (vtk_data_type)
395 {
396 case VTK_UNIFORM_GRID:
397 dataObject = this->RequestImageData(xmfGrid, true);
398 break;
399
400 case VTK_IMAGE_DATA:
401 dataObject = this->RequestImageData(xmfGrid, false);
402 break;
403
404 case VTK_STRUCTURED_GRID:
405 dataObject = this->RequestStructuredGrid(xmfGrid);
406 break;
407
408 case VTK_RECTILINEAR_GRID:
409 dataObject = this->RequestRectilinearGrid(xmfGrid);
410 break;
411
412 case VTK_UNSTRUCTURED_GRID:
413 dataObject = this->ReadUnstructuredGrid(xmfGrid);
414 break;
415
416 default:
417 // un-handled case.
418 return nullptr;
419 }
420
421 if (caching)
422 {
423 cache[blockId].dataset = vtkDataSet::SafeDownCast(dataObject);
424 dataObject->Register(nullptr);
425 }
426 return dataObject;
427 }
428
429 //------------------------------------------------------------------------------
GetNumberOfPointsPerCell(int vtk_cell_type)430 int vtkXdmfHeavyData::GetNumberOfPointsPerCell(int vtk_cell_type)
431 {
432 switch (vtk_cell_type)
433 {
434 case VTK_POLY_VERTEX:
435 return 0;
436 case VTK_POLY_LINE:
437 return 0;
438 case VTK_POLYGON:
439 return 0;
440
441 case VTK_TRIANGLE:
442 return 3;
443 case VTK_QUAD:
444 return 4;
445 case VTK_TETRA:
446 return 4;
447 case VTK_PYRAMID:
448 return 5;
449 case VTK_WEDGE:
450 return 6;
451 case VTK_HEXAHEDRON:
452 return 8;
453 case VTK_QUADRATIC_EDGE:
454 return 3;
455 case VTK_QUADRATIC_TRIANGLE:
456 return 6;
457 case VTK_QUADRATIC_QUAD:
458 return 8;
459 case VTK_BIQUADRATIC_QUAD:
460 return 9;
461 case VTK_QUADRATIC_TETRA:
462 return 10;
463 case VTK_QUADRATIC_PYRAMID:
464 return 13;
465 case VTK_QUADRATIC_WEDGE:
466 return 15;
467 case VTK_BIQUADRATIC_QUADRATIC_WEDGE:
468 return 18;
469 case VTK_QUADRATIC_HEXAHEDRON:
470 return 20;
471 case VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON:
472 return 24;
473 case VTK_TRIQUADRATIC_HEXAHEDRON:
474 return 24;
475 }
476 return -1;
477 }
478 //------------------------------------------------------------------------------
GetVTKCellType(XdmfInt32 topologyType)479 int vtkXdmfHeavyData::GetVTKCellType(XdmfInt32 topologyType)
480 {
481 switch (topologyType)
482 {
483 case XDMF_POLYVERTEX:
484 return VTK_POLY_VERTEX;
485 case XDMF_POLYLINE:
486 return VTK_POLY_LINE;
487 case XDMF_POLYGON:
488 return VTK_POLYGON; // FIXME: should this not be treated as mixed?
489 case XDMF_TRI:
490 return VTK_TRIANGLE;
491 case XDMF_QUAD:
492 return VTK_QUAD;
493 case XDMF_TET:
494 return VTK_TETRA;
495 case XDMF_PYRAMID:
496 return VTK_PYRAMID;
497 case XDMF_WEDGE:
498 return VTK_WEDGE;
499 case XDMF_HEX:
500 return VTK_HEXAHEDRON;
501 case XDMF_EDGE_3:
502 return VTK_QUADRATIC_EDGE;
503 case XDMF_TRI_6:
504 return VTK_QUADRATIC_TRIANGLE;
505 case XDMF_QUAD_8:
506 return VTK_QUADRATIC_QUAD;
507 case XDMF_QUAD_9:
508 return VTK_BIQUADRATIC_QUAD;
509 case XDMF_TET_10:
510 return VTK_QUADRATIC_TETRA;
511 case XDMF_PYRAMID_13:
512 return VTK_QUADRATIC_PYRAMID;
513 case XDMF_WEDGE_15:
514 return VTK_QUADRATIC_WEDGE;
515 case XDMF_WEDGE_18:
516 return VTK_BIQUADRATIC_QUADRATIC_WEDGE;
517 case XDMF_HEX_20:
518 return VTK_QUADRATIC_HEXAHEDRON;
519 case XDMF_HEX_24:
520 return VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON;
521 case XDMF_HEX_27:
522 return VTK_TRIQUADRATIC_HEXAHEDRON;
523 case XDMF_MIXED:
524 return VTK_NUMBER_OF_CELL_TYPES;
525 }
526 // XdmfErrorMessage("Unknown Topology Type = "
527 // << xmfGrid->GetTopology()->GetTopologyType());
528 return VTK_EMPTY_CELL;
529 }
530
531 //------------------------------------------------------------------------------
ReadUnstructuredGrid(XdmfGrid * xmfGrid)532 vtkDataObject* vtkXdmfHeavyData::ReadUnstructuredGrid(XdmfGrid* xmfGrid)
533 {
534 vtkSmartPointer<vtkUnstructuredGrid> ugData = vtkSmartPointer<vtkUnstructuredGrid>::New();
535
536 // BUG #12527. For non-partitioned data, don't read unstructured grid on
537 // process id > 0.
538 if (this->Piece != 0 && this->Domain->GetNumberOfGrids() == 1 &&
539 this->Domain->GetVTKDataType() == VTK_UNSTRUCTURED_GRID &&
540 this->Domain->GetSetsSelection()->GetNumberOfArrays() == 0)
541 {
542 ugData->Register(nullptr);
543 return ugData;
544 }
545
546 XdmfTopology* xmfTopology = xmfGrid->GetTopology();
547 XdmfArray* xmfConnectivity = xmfTopology->GetConnectivity();
548
549 int vtk_cell_type = vtkXdmfHeavyData::GetVTKCellType(xmfTopology->GetTopologyType());
550
551 if (vtk_cell_type == VTK_EMPTY_CELL)
552 {
553 // invalid topology.
554 return nullptr;
555 }
556
557 if (vtk_cell_type != VTK_NUMBER_OF_CELL_TYPES)
558 // i.e. topologyType != XDMF_MIXED
559 {
560 // all cells are of the same type.
561 XdmfInt32 numPointsPerCell = xmfTopology->GetNodesPerElement();
562
563 // FIXME: is this needed, shouldn't xmfTopology->GetNodesPerElement()
564 // return the correct value always?
565 if (xmfConnectivity->GetRank() == 2)
566 {
567 numPointsPerCell = xmfConnectivity->GetDimension(1);
568 }
569
570 /* Create Cell Type Array */
571 XdmfInt64 conn_length = xmfConnectivity->GetNumberOfElements();
572 std::vector<XdmfInt64> xmfConnections(conn_length);
573 xmfConnectivity->GetValues(0, xmfConnections.data(), conn_length);
574
575 vtkIdType numCells = xmfTopology->GetShapeDesc()->GetNumberOfElements();
576
577 vtkNew<vtkIdTypeArray> conn;
578 vtkNew<vtkIdTypeArray> offsets;
579
580 offsets->SetNumberOfTuples(numCells + 1);
581
582 { // Fill offsets: {0, 1 * cellSize, 2 * cellSize, ..., numCells * cellSize}
583 vtkIdType offset = -numPointsPerCell;
584 auto generator = [&]() -> vtkIdType { return offset += numPointsPerCell; };
585 auto range = vtk::DataArrayValueRange<1>(offsets);
586 std::generate(range.begin(), range.end(), generator);
587 }
588
589 conn->SetNumberOfTuples(numPointsPerCell * numCells);
590
591 // Fill connections (just copy xmfConnections)
592 { // Need to convert explicitly to silence warnings:
593 auto range = vtk::DataArrayValueRange<1>(conn);
594 std::transform(xmfConnections.cbegin(),
595 xmfConnections.cbegin() + (numPointsPerCell * numCells), range.begin(),
596 [](XdmfInt64 val) -> vtkIdType { return static_cast<vtkIdType>(val); });
597 }
598
599 // Construct and set the cell array
600 vtkNew<vtkCellArray> cells;
601 cells->SetData(offsets, conn);
602 ugData->SetCells(vtk_cell_type, cells);
603 }
604 else
605 {
606 // We have cells with mixed types.
607 XdmfInt64 conn_length = xmfGrid->GetTopology()->GetConnectivity()->GetNumberOfElements();
608 std::vector<XdmfInt64> xmfConnections(static_cast<size_t>(conn_length));
609 xmfConnectivity->GetValues(0, xmfConnections.data(), conn_length);
610
611 vtkIdType numCells = xmfTopology->GetShapeDesc()->GetNumberOfElements();
612 vtkNew<vtkUnsignedCharArray> cell_types;
613 cell_types->SetNumberOfTuples(numCells);
614
615 vtkNew<vtkIdTypeArray> offsets;
616 offsets->SetNumberOfTuples(numCells + 1);
617
618 vtkNew<vtkIdTypeArray> conn;
619 // This may be an overestimate; will correct after filling.
620 conn->SetNumberOfTuples(static_cast<vtkIdType>(conn_length));
621
622 vtkIdType offset = 0;
623 vtkIdType index = 0;
624 vtkIdType connIndex = 0;
625 for (vtkIdType cc = 0; cc < numCells; cc++)
626 {
627 int vtk_cell_typeI = this->GetVTKCellType(xmfConnections[index++]);
628 XdmfInt32 numPointsPerCell = this->GetNumberOfPointsPerCell(vtk_cell_typeI);
629 if (numPointsPerCell == -1)
630 {
631 // encountered an unknown cell.
632 return nullptr;
633 }
634
635 if (numPointsPerCell == 0)
636 {
637 // cell type does not have a fixed number of points in which case the
638 // next entry in xmfConnections tells us the number of points.
639 numPointsPerCell = xmfConnections[index++];
640 }
641
642 cell_types->SetValue(cc, static_cast<unsigned char>(vtk_cell_typeI));
643 offsets->SetValue(cc, offset);
644 offset += numPointsPerCell;
645
646 for (vtkIdType i = 0; i < numPointsPerCell; i++)
647 {
648 conn->SetValue(connIndex++, xmfConnections[index++]);
649 }
650 }
651 offsets->SetValue(numCells, offset); // final offset value
652
653 // Resize the Array to the Proper Size
654 conn->Resize(connIndex);
655
656 // Create and set the cell array:
657 vtkNew<vtkCellArray> cells;
658 cells->SetData(offsets, conn);
659 ugData->SetCells(cell_types, cells);
660 }
661
662 // Read the geometry.
663 vtkPoints* points = this->ReadPoints(xmfGrid->GetGeometry());
664 if (!points)
665 {
666 // failed to read points.
667 return nullptr;
668 }
669 ugData->SetPoints(points);
670 points->Delete();
671
672 this->ReadAttributes(ugData, xmfGrid);
673
674 // Read ghost cell/point information.
675 this->ReadGhostSets(ugData, xmfGrid);
676
677 // If this grid has sets defined on it, then we need to read those as well
678 vtkMultiBlockDataSet* sets = this->ReadSets(ugData, xmfGrid);
679 if (sets)
680 {
681 return sets;
682 }
683
684 ugData->Register(nullptr);
685 return ugData;
686 }
687
vtkExtentsAreValid(int exts[6])688 inline bool vtkExtentsAreValid(int exts[6])
689 {
690 return exts[1] >= exts[0] && exts[3] >= exts[2] && exts[5] >= exts[4];
691 }
692
vtkExtentsAreEqual(int * exts1,int * exts2)693 inline bool vtkExtentsAreEqual(int* exts1, int* exts2)
694 {
695 if (!exts1 && !exts2)
696 {
697 return true;
698 }
699 if (!exts1 || !exts2)
700 {
701 return false;
702 }
703 return (exts1[0] == exts2[0] && exts1[1] == exts2[1] && exts1[2] == exts2[2] &&
704 exts1[3] == exts2[3] && exts1[4] == exts2[4] && exts1[5] == exts2[5]);
705 }
706
707 //------------------------------------------------------------------------------
RequestRectilinearGrid(XdmfGrid * xmfGrid)708 vtkRectilinearGrid* vtkXdmfHeavyData::RequestRectilinearGrid(XdmfGrid* xmfGrid)
709 {
710 vtkSmartPointer<vtkRectilinearGrid> rg = vtkSmartPointer<vtkRectilinearGrid>::New();
711 int whole_extents[6];
712 int update_extents[6];
713 this->Domain->GetWholeExtent(xmfGrid, whole_extents);
714
715 if (!vtkExtentsAreValid(this->Extents))
716 {
717 // if this->Extents are not valid, then simply read the whole image.
718 memcpy(update_extents, whole_extents, sizeof(int) * 6);
719 }
720 else
721 {
722 memcpy(update_extents, this->Extents, sizeof(int) * 6);
723 }
724
725 // convert to stridden update extents.
726 int scaled_extents[6];
727 vtkScaleExtents(update_extents, scaled_extents, this->Stride);
728 int scaled_dims[3];
729 vtkGetDims(scaled_extents, scaled_dims);
730
731 rg->SetExtent(scaled_extents);
732
733 // Now read rectilinear geometry.
734 XdmfGeometry* xmfGeometry = xmfGrid->GetGeometry();
735
736 vtkSmartPointer<vtkDoubleArray> xarray = vtkSmartPointer<vtkDoubleArray>::New();
737 xarray->SetNumberOfTuples(scaled_dims[0]);
738
739 vtkSmartPointer<vtkDoubleArray> yarray = vtkSmartPointer<vtkDoubleArray>::New();
740 yarray->SetNumberOfTuples(scaled_dims[1]);
741
742 vtkSmartPointer<vtkDoubleArray> zarray = vtkSmartPointer<vtkDoubleArray>::New();
743 zarray->SetNumberOfTuples(scaled_dims[2]);
744
745 rg->SetXCoordinates(xarray);
746 rg->SetYCoordinates(yarray);
747 rg->SetZCoordinates(zarray);
748
749 switch (xmfGeometry->GetGeometryType())
750 {
751 case XDMF_GEOMETRY_ORIGIN_DXDY:
752 case XDMF_GEOMETRY_ORIGIN_DXDYDZ:
753 {
754 XdmfFloat64* origin = xmfGeometry->GetOrigin();
755 XdmfFloat64* dxdydz = xmfGeometry->GetDxDyDz();
756 for (int cc = scaled_extents[0]; cc <= scaled_extents[1]; cc++)
757 {
758 xarray->GetPointer(0)[cc - scaled_extents[0]] =
759 origin[0] + (dxdydz[0] * cc * this->Stride[0]);
760 }
761 for (int cc = scaled_extents[2]; cc <= scaled_extents[3]; cc++)
762 {
763 yarray->GetPointer(0)[cc - scaled_extents[2]] =
764 origin[1] + (dxdydz[1] * cc * this->Stride[1]);
765 }
766 for (int cc = scaled_extents[4]; cc <= scaled_extents[5]; cc++)
767 {
768 zarray->GetPointer(0)[cc - scaled_extents[4]] =
769 origin[2] + (dxdydz[2] * cc * this->Stride[2]);
770 }
771 }
772 break;
773
774 case XDMF_GEOMETRY_VXVY:
775 {
776 // Note:
777 // XDMF and VTK structured extents are reversed
778 // Where I varies fastest, VTK's convention is IJK, but XDMF's is KJI
779 // However, users naturally don't want VXVY to mean VZVY.
780 // Let's accept VisIt's interpretation of this 2D case
781 // (KJI is ZXY where Z=0).
782 xarray->SetNumberOfTuples(scaled_dims[1]);
783 yarray->SetNumberOfTuples(scaled_dims[2]);
784 zarray->SetNumberOfTuples(scaled_dims[0]);
785 rg->SetExtent(scaled_extents[2], scaled_extents[3], scaled_extents[4], scaled_extents[5],
786 scaled_extents[0], scaled_extents[1]);
787 xmfGeometry->GetVectorX()->GetValues(
788 update_extents[2], xarray->GetPointer(0), scaled_dims[1], this->Stride[1]);
789 xmfGeometry->GetVectorY()->GetValues(
790 update_extents[4], yarray->GetPointer(0), scaled_dims[2], this->Stride[2]);
791 zarray->FillComponent(0, 0);
792 }
793 break;
794
795 case XDMF_GEOMETRY_VXVYVZ:
796 {
797 xmfGeometry->GetVectorX()->GetValues(
798 update_extents[0], xarray->GetPointer(0), scaled_dims[0], this->Stride[0]);
799 xmfGeometry->GetVectorY()->GetValues(
800 update_extents[2], yarray->GetPointer(0), scaled_dims[1], this->Stride[1]);
801 xmfGeometry->GetVectorZ()->GetValues(
802 update_extents[4], zarray->GetPointer(0), scaled_dims[2], this->Stride[2]);
803 }
804 break;
805
806 default:
807 vtkErrorWithObjectMacro(this->Reader,
808 "Geometry type : " << xmfGeometry->GetGeometryTypeAsString() << " is not supported for "
809 << xmfGrid->GetTopology()->GetTopologyTypeAsString());
810 return nullptr;
811 }
812
813 this->ReadAttributes(rg, xmfGrid, update_extents);
814 rg->Register(nullptr);
815 return rg;
816 }
817
818 //------------------------------------------------------------------------------
RequestStructuredGrid(XdmfGrid * xmfGrid)819 vtkStructuredGrid* vtkXdmfHeavyData::RequestStructuredGrid(XdmfGrid* xmfGrid)
820 {
821 vtkStructuredGrid* sg = vtkStructuredGrid::New();
822
823 int whole_extents[6];
824 int update_extents[6];
825 this->Domain->GetWholeExtent(xmfGrid, whole_extents);
826
827 if (!vtkExtentsAreValid(this->Extents))
828 {
829 // if this->Extents are not valid, then simply read the whole image.
830 memcpy(update_extents, whole_extents, sizeof(int) * 6);
831 }
832 else
833 {
834 memcpy(update_extents, this->Extents, sizeof(int) * 6);
835 }
836
837 int scaled_extents[6];
838 vtkScaleExtents(update_extents, scaled_extents, this->Stride);
839 sg->SetExtent(scaled_extents);
840
841 vtkPoints* points = this->ReadPoints(xmfGrid->GetGeometry(), update_extents, whole_extents);
842 sg->SetPoints(points);
843 points->Delete();
844
845 this->ReadAttributes(sg, xmfGrid, update_extents);
846 return sg;
847 }
848
849 //------------------------------------------------------------------------------
RequestImageData(XdmfGrid * xmfGrid,bool use_uniform_grid)850 vtkImageData* vtkXdmfHeavyData::RequestImageData(XdmfGrid* xmfGrid, bool use_uniform_grid)
851 {
852 vtkImageData* imageData =
853 use_uniform_grid ? static_cast<vtkImageData*>(vtkUniformGrid::New()) : vtkImageData::New();
854
855 int whole_extents[6];
856 this->Domain->GetWholeExtent(xmfGrid, whole_extents);
857
858 int update_extents[6];
859
860 if (!vtkExtentsAreValid(this->Extents))
861 {
862 // if this->Extents are not valid, then simply read the whole image.
863 memcpy(update_extents, whole_extents, sizeof(int) * 6);
864 }
865 else
866 {
867 memcpy(update_extents, this->Extents, sizeof(int) * 6);
868 }
869
870 int scaled_extents[6];
871 vtkScaleExtents(update_extents, scaled_extents, this->Stride);
872 imageData->SetExtent(scaled_extents);
873
874 double origin[3], spacing[3];
875 if (!this->Domain->GetOriginAndSpacing(xmfGrid, origin, spacing))
876 {
877 vtkErrorWithObjectMacro(this->Reader,
878 "Could not determine image-data origin and spacing. "
879 "Required geometry type is ORIGIN_DXDY or ORIGIN_DXDYDZ. "
880 "The specified geometry type is : "
881 << xmfGrid->GetGeometry()->GetGeometryTypeAsString());
882 // release image data.
883 imageData->Delete();
884 return nullptr;
885 }
886 imageData->SetOrigin(origin);
887 imageData->SetSpacing(
888 spacing[0] * this->Stride[0], spacing[1] * this->Stride[1], spacing[2] * this->Stride[2]);
889 this->ReadAttributes(imageData, xmfGrid, update_extents);
890 return imageData;
891 }
892
893 //------------------------------------------------------------------------------
ReadPoints(XdmfGeometry * xmfGeometry,int * update_extents,int * whole_extents)894 vtkPoints* vtkXdmfHeavyData::ReadPoints(
895 XdmfGeometry* xmfGeometry, int* update_extents /*=nullptr*/, int* whole_extents /*=nullptr*/)
896 {
897 XdmfInt32 geomType = xmfGeometry->GetGeometryType();
898
899 if (geomType != XDMF_GEOMETRY_X_Y_Z && geomType != XDMF_GEOMETRY_XYZ &&
900 geomType != XDMF_GEOMETRY_X_Y && geomType != XDMF_GEOMETRY_XY)
901 {
902 return nullptr;
903 }
904
905 XdmfArray* xmfPoints = xmfGeometry->GetPoints();
906 if (!xmfPoints)
907 {
908 XdmfErrorMessage("No Points to Set");
909 return nullptr;
910 }
911
912 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
913
914 if (xmfPoints->GetNumberType() == XDMF_FLOAT32_TYPE)
915 {
916 vtkFloatArray* da = vtkFloatArray::New();
917 da->SetNumberOfComponents(3);
918 points->SetData(da);
919 da->Delete();
920 }
921 else // means == XDMF_FLOAT64_TYPE
922 {
923 vtkDoubleArray* da = vtkDoubleArray::New();
924 da->SetNumberOfComponents(3);
925 points->SetData(da);
926 da->Delete();
927 }
928
929 XdmfInt64 numGeometryPoints = xmfGeometry->GetNumberOfPoints();
930 vtkIdType numPoints = numGeometryPoints;
931 bool structured_data = false;
932 if (update_extents && whole_extents)
933 {
934 // we are reading a sub-extent.
935 structured_data = true;
936 int scaled_extents[6];
937 int scaled_dims[3];
938 vtkScaleExtents(update_extents, scaled_extents, this->Stride);
939 vtkGetDims(scaled_extents, scaled_dims);
940 numPoints = (scaled_dims[0] * scaled_dims[1] * scaled_dims[2]);
941 }
942 points->SetNumberOfPoints(numPoints);
943
944 if (!structured_data)
945 {
946 // read all the points.
947 switch (points->GetData()->GetDataType())
948 {
949 case VTK_DOUBLE:
950 xmfPoints->GetValues(
951 0, reinterpret_cast<double*>(points->GetVoidPointer(0)), numPoints * 3);
952 break;
953
954 case VTK_FLOAT:
955 xmfPoints->GetValues(0, reinterpret_cast<float*>(points->GetVoidPointer(0)), numPoints * 3);
956 break;
957
958 default:
959 return nullptr;
960 }
961 }
962 else
963 {
964 // treating the points as structured points
965 std::vector<XdmfFloat64> tempPoints(numGeometryPoints * 3);
966 xmfPoints->GetValues(0, tempPoints.data(), numGeometryPoints * 3);
967 vtkIdType pointId = 0;
968 int xdmf_dims[3];
969 vtkGetDims(whole_extents, xdmf_dims);
970
971 for (int z = update_extents[4]; z <= update_extents[5]; z++)
972 {
973 if ((z - update_extents[4]) % this->Stride[2])
974 {
975 continue;
976 }
977
978 for (int y = update_extents[2]; y <= update_extents[3]; y++)
979 {
980 if ((y - update_extents[2]) % this->Stride[1])
981 {
982 continue;
983 }
984
985 for (int x = update_extents[0]; x <= update_extents[1]; x++)
986 {
987 if ((x - update_extents[0]) % this->Stride[0])
988 {
989 continue;
990 }
991
992 int xdmf_index[3] = { x, y, z };
993 XdmfInt64 offset = vtkStructuredData::ComputePointId(xdmf_dims, xdmf_index);
994 points->SetPoint(pointId, tempPoints[3 * offset], tempPoints[3 * offset + 1],
995 tempPoints[3 * offset + 2]);
996 pointId++;
997 }
998 }
999 }
1000 }
1001
1002 points->Register(nullptr);
1003 return points;
1004 }
1005
1006 //------------------------------------------------------------------------------
ReadAttributes(vtkDataSet * dataSet,XdmfGrid * xmfGrid,int * update_extents)1007 bool vtkXdmfHeavyData::ReadAttributes(vtkDataSet* dataSet, XdmfGrid* xmfGrid, int* update_extents)
1008 {
1009 int data_dimensionality = vtkXdmfDomain::GetDataDimensionality(xmfGrid);
1010
1011 int numAttributes = xmfGrid->GetNumberOfAttributes();
1012 for (int cc = 0; cc < numAttributes; cc++)
1013 {
1014 XdmfAttribute* xmfAttribute = xmfGrid->GetAttribute(cc);
1015 const char* attrName = xmfAttribute->GetName();
1016 int attrCenter = xmfAttribute->GetAttributeCenter();
1017 if (!attrName)
1018 {
1019 vtkWarningWithObjectMacro(this->Reader, "Skipping unnamed attributes.");
1020 continue;
1021 }
1022
1023 vtkFieldData* fieldData = nullptr;
1024 // skip disabled arrays.
1025 switch (attrCenter)
1026 {
1027 case XDMF_ATTRIBUTE_CENTER_GRID:
1028 fieldData = dataSet->GetFieldData();
1029 break;
1030
1031 case XDMF_ATTRIBUTE_CENTER_CELL:
1032 if (!this->Domain->GetCellArraySelection()->ArrayIsEnabled(attrName))
1033 {
1034 continue;
1035 }
1036 fieldData = dataSet->GetCellData();
1037 break;
1038
1039 case XDMF_ATTRIBUTE_CENTER_NODE:
1040 if (!this->Domain->GetPointArraySelection()->ArrayIsEnabled(attrName))
1041 {
1042 continue;
1043 }
1044 fieldData = dataSet->GetPointData();
1045 break;
1046
1047 case XDMF_ATTRIBUTE_CENTER_FACE:
1048 case XDMF_ATTRIBUTE_CENTER_EDGE:
1049 default:
1050 vtkWarningWithObjectMacro(this->Reader,
1051 "Skipping attribute " << attrName << " at "
1052 << xmfAttribute->GetAttributeCenterAsString());
1053 continue; // unhandled.
1054 }
1055
1056 vtkDataArray* array = this->ReadAttribute(xmfAttribute, data_dimensionality, update_extents);
1057 if (array)
1058 {
1059 array->SetName(attrName);
1060 fieldData->AddArray(array);
1061 bool is_active = xmfAttribute->GetActive() != 0;
1062 vtkDataSetAttributes* attributes = vtkDataSetAttributes::SafeDownCast(fieldData);
1063 if (attributes)
1064 {
1065 // make attribute active.
1066 switch (xmfAttribute->GetAttributeType())
1067 {
1068 case XDMF_ATTRIBUTE_TYPE_SCALAR:
1069 if (is_active || attributes->GetScalars() == nullptr)
1070 {
1071 attributes->SetActiveScalars(attrName);
1072 }
1073 break;
1074
1075 case XDMF_ATTRIBUTE_TYPE_VECTOR:
1076 if (is_active || attributes->GetVectors() == nullptr)
1077 {
1078 attributes->SetActiveVectors(attrName);
1079 }
1080 break;
1081
1082 case XDMF_ATTRIBUTE_TYPE_TENSOR:
1083 case XDMF_ATTRIBUTE_TYPE_TENSOR6:
1084 if (is_active || attributes->GetTensors() == nullptr)
1085 {
1086 attributes->SetActiveTensors(attrName);
1087 }
1088 break;
1089
1090 case XDMF_ATTRIBUTE_TYPE_GLOBALID:
1091 if (is_active || attributes->GetGlobalIds() == nullptr)
1092 {
1093 attributes->SetActiveGlobalIds(attrName);
1094 }
1095 }
1096 }
1097 array->Delete();
1098 }
1099 }
1100 return true;
1101 }
1102
1103 // used to convert a symmetric tensor to a regular tensor.
1104 template <class T>
vtkConvertTensor6(T * source,T * dest,vtkIdType numTensors)1105 void vtkConvertTensor6(T* source, T* dest, vtkIdType numTensors)
1106 {
1107 for (vtkIdType cc = 0; cc < numTensors; cc++)
1108 {
1109 dest[cc * 9 + 0] = source[cc * 6 + 0];
1110 dest[cc * 9 + 1] = source[cc * 6 + 1];
1111 dest[cc * 9 + 2] = source[cc * 6 + 2];
1112
1113 dest[cc * 9 + 3] = source[cc * 6 + 1];
1114 dest[cc * 9 + 4] = source[cc * 6 + 3];
1115 dest[cc * 9 + 5] = source[cc * 6 + 4];
1116
1117 dest[cc * 9 + 6] = source[cc * 6 + 2];
1118 dest[cc * 9 + 7] = source[cc * 6 + 4];
1119 dest[cc * 9 + 8] = source[cc * 6 + 5];
1120 }
1121 }
1122
1123 //------------------------------------------------------------------------------
ReadAttribute(XdmfAttribute * xmfAttribute,int data_dimensionality,int * update_extents)1124 vtkDataArray* vtkXdmfHeavyData::ReadAttribute(
1125 XdmfAttribute* xmfAttribute, int data_dimensionality, int* update_extents /*=0*/)
1126 {
1127 if (!xmfAttribute)
1128 {
1129 return nullptr;
1130 }
1131
1132 int attrType = xmfAttribute->GetAttributeType();
1133 int attrCenter = xmfAttribute->GetAttributeCenter();
1134 int numComponents = 1;
1135
1136 XdmfDataItem xmfDataItem;
1137 xmfDataItem.SetDOM(xmfAttribute->GetDOM());
1138 xmfDataItem.SetElement(xmfAttribute->GetDOM()->FindDataElement(0, xmfAttribute->GetElement()));
1139 xmfDataItem.UpdateInformation();
1140
1141 XdmfInt64 data_dims[XDMF_MAX_DIMENSION];
1142 int data_rank = xmfDataItem.GetDataDesc()->GetShape(data_dims);
1143
1144 switch (attrType)
1145 {
1146 case XDMF_ATTRIBUTE_TYPE_TENSOR:
1147 numComponents = 9;
1148 break;
1149 case XDMF_ATTRIBUTE_TYPE_TENSOR6:
1150 numComponents = 6;
1151 break;
1152 case XDMF_ATTRIBUTE_TYPE_VECTOR:
1153 numComponents = 3;
1154 break;
1155 case XDMF_ATTRIBUTE_TYPE_MATRIX:
1156 numComponents = data_dims[data_rank - 1];
1157 break;
1158 default:
1159 numComponents = 1;
1160 break;
1161 }
1162
1163 // Handle 2D vectors
1164 if (attrType == XDMF_ATTRIBUTE_TYPE_VECTOR && data_dims[data_rank - 1] == 2)
1165 {
1166 numComponents = 2;
1167 }
1168
1169 if (update_extents && attrCenter != XDMF_ATTRIBUTE_CENTER_GRID)
1170 {
1171 // for hyperslab selection to work, the data shape must match the topology
1172 // shape.
1173 if (data_rank < 0)
1174 {
1175 vtkErrorWithObjectMacro(this->Reader, "Unsupported attribute rank: " << data_rank);
1176 return nullptr;
1177 }
1178 if (data_rank > (data_dimensionality + 1))
1179 {
1180 vtkErrorWithObjectMacro(
1181 this->Reader, "The data_dimensionality and topology dimensionality mismatch");
1182 return nullptr;
1183 }
1184 XdmfInt64 start[4] = { update_extents[4], update_extents[2], update_extents[0], 0 };
1185 XdmfInt64 stride[4] = { this->Stride[2], this->Stride[1], this->Stride[0], 1 };
1186 XdmfInt64 count[4] = { 0, 0, 0, 0 };
1187 int scaled_dims[3];
1188 int scaled_extents[6];
1189 vtkScaleExtents(update_extents, scaled_extents, this->Stride);
1190 vtkGetDims(scaled_extents, scaled_dims);
1191 count[0] = (scaled_dims[2] - 1);
1192 count[1] = (scaled_dims[1] - 1);
1193 count[2] = (scaled_dims[0] - 1);
1194 if (data_rank == (data_dimensionality + 1))
1195 {
1196 // this refers the number of components in the attribute.
1197 count[data_dimensionality] = data_dims[data_dimensionality];
1198 }
1199
1200 if (attrCenter == XDMF_ATTRIBUTE_CENTER_NODE)
1201 {
1202 // Point count is 1 + cell extent if not a single layer
1203 count[0] += 1; //((update_extents[5] - update_extents[4]) > 0)? 1 : 0;
1204 count[1] += 1; //((update_extents[3] - update_extents[2]) > 0)? 1 : 0;
1205 count[2] += 1; //((update_extents[1] - update_extents[0]) > 0)? 1 : 0;
1206 }
1207 xmfDataItem.GetDataDesc()->SelectHyperSlab(start, stride, count);
1208 }
1209
1210 if (xmfDataItem.Update() == XDMF_FAIL)
1211 {
1212 vtkErrorWithObjectMacro(this->Reader, "Failed to read attribute data");
1213 return nullptr;
1214 }
1215
1216 vtkXdmfDataArray* xmfConvertor = vtkXdmfDataArray::New();
1217 vtkDataArray* dataArray = xmfConvertor->FromXdmfArray(
1218 xmfDataItem.GetArray()->GetTagName(), 1, data_rank, numComponents, 0);
1219 xmfConvertor->Delete();
1220
1221 if (attrType == XDMF_ATTRIBUTE_TYPE_TENSOR6)
1222 {
1223 // convert Tensor6 to Tensor.
1224 vtkDataArray* tensor = dataArray->NewInstance();
1225 vtkIdType numTensors = dataArray->GetNumberOfTuples();
1226 tensor->SetNumberOfComponents(9);
1227 tensor->SetNumberOfTuples(numTensors);
1228
1229 // Copy Symmetrical Tensor Values to Correct Positions in 3x3 matrix
1230 void* source = dataArray->GetVoidPointer(0);
1231 void* dest = tensor->GetVoidPointer(0);
1232 switch (tensor->GetDataType())
1233 {
1234 vtkTemplateMacro(vtkConvertTensor6(
1235 reinterpret_cast<VTK_TT*>(source), reinterpret_cast<VTK_TT*>(dest), numTensors));
1236 }
1237 dataArray->Delete();
1238 return tensor;
1239 }
1240
1241 if (attrType == XDMF_ATTRIBUTE_TYPE_VECTOR && numComponents == 2)
1242 {
1243 // convert 2D vectors to 3-tuple vectors with 0.0 in the z component
1244 vtkDataArray* vector3D = dataArray->NewInstance();
1245 vtkIdType numVectors = dataArray->GetNumberOfTuples();
1246 vector3D->SetNumberOfComponents(3);
1247 vector3D->SetNumberOfTuples(numVectors);
1248
1249 // Add 0.0 to third component of vector
1250 auto inputRange = vtk::DataArrayTupleRange<2>(dataArray);
1251 auto outputRange = vtk::DataArrayTupleRange<3>(vector3D);
1252 for (auto i = 0; i < inputRange.size(); ++i)
1253 {
1254 outputRange[i][0] = inputRange[i][0];
1255 outputRange[i][1] = inputRange[i][1];
1256 outputRange[i][2] = 0.0;
1257 }
1258 dataArray->Delete();
1259 return vector3D;
1260 }
1261
1262 return dataArray;
1263 }
1264
1265 //------------------------------------------------------------------------------
1266 // Read ghost cell/point information. This is simply loaded info a
1267 // vtkGhostType attribute array.
ReadGhostSets(vtkDataSet * dataSet,XdmfGrid * xmfGrid,int * vtkNotUsed (update_extents))1268 bool vtkXdmfHeavyData::ReadGhostSets(
1269 vtkDataSet* dataSet, XdmfGrid* xmfGrid, int* vtkNotUsed(update_extents) /*=0*/)
1270 {
1271 // int data_dimensionality = this->Domain->GetDataDimensionality(xmfGrid);
1272 for (int cc = 0; cc < xmfGrid->GetNumberOfSets(); cc++)
1273 {
1274 XdmfSet* xmfSet = xmfGrid->GetSets(cc);
1275 int ghost_value = xmfSet->GetGhost();
1276 if (ghost_value <= 0)
1277 {
1278 // not a ghost-set, simply continue.
1279 continue;
1280 }
1281 XdmfInt32 setCenter = xmfSet->GetSetType();
1282 vtkIdType numElems = 0;
1283 vtkDataSetAttributes* dsa = nullptr;
1284 unsigned char ghostFlag = 0;
1285 switch (setCenter)
1286 {
1287 case XDMF_SET_TYPE_NODE:
1288 dsa = dataSet->GetPointData();
1289 numElems = dataSet->GetNumberOfPoints();
1290 ghostFlag = vtkDataSetAttributes::DUPLICATEPOINT;
1291 break;
1292
1293 case XDMF_SET_TYPE_CELL:
1294 dsa = dataSet->GetCellData();
1295 numElems = dataSet->GetNumberOfCells();
1296 ghostFlag = vtkDataSetAttributes::DUPLICATECELL;
1297 break;
1298
1299 default:
1300 vtkWarningWithObjectMacro(
1301 this->Reader, "Only ghost-cells and ghost-nodes are currently supported.");
1302 continue;
1303 }
1304
1305 vtkUnsignedCharArray* ghosts =
1306 vtkArrayDownCast<vtkUnsignedCharArray>(dsa->GetArray(vtkDataSetAttributes::GhostArrayName()));
1307 if (!ghosts)
1308 {
1309 ghosts = vtkUnsignedCharArray::New();
1310 ghosts->SetName(vtkDataSetAttributes::GhostArrayName());
1311 ghosts->SetNumberOfComponents(1);
1312 ghosts->SetNumberOfTuples(numElems);
1313 ghosts->FillComponent(0, 0);
1314 dsa->AddArray(ghosts);
1315 ghosts->Delete();
1316 }
1317
1318 unsigned char* ptrGhosts = ghosts->GetPointer(0);
1319
1320 // Read heavy data. We cannot do anything smart if update_extents or stride
1321 // is specified here. We have to read the entire set and then prune it.
1322 xmfSet->Update();
1323
1324 XdmfArray* xmfIds = xmfSet->GetIds();
1325 XdmfInt64 numIds = xmfIds->GetNumberOfElements();
1326 std::vector<XdmfInt64> ids(numIds + 1);
1327 xmfIds->GetValues(0, ids.data(), numIds);
1328
1329 // release the heavy data that was read.
1330 xmfSet->Release();
1331
1332 for (vtkIdType kk = 0; kk < numIds; kk++)
1333 {
1334 if (ids[kk] < 0 || ids[kk] > numElems)
1335 {
1336 vtkWarningWithObjectMacro(this->Reader, "No such cell or point exists: " << ids[kk]);
1337 continue;
1338 }
1339 ptrGhosts[ids[kk]] = ghostFlag;
1340 }
1341 }
1342 return true;
1343 }
1344
1345 //------------------------------------------------------------------------------
ReadSets(vtkDataSet * dataSet,XdmfGrid * xmfGrid,int * vtkNotUsed (update_extents))1346 vtkMultiBlockDataSet* vtkXdmfHeavyData::ReadSets(
1347 vtkDataSet* dataSet, XdmfGrid* xmfGrid, int* vtkNotUsed(update_extents) /*=0*/)
1348 {
1349 unsigned int number_of_sets = 0;
1350 for (int cc = 0; cc < xmfGrid->GetNumberOfSets(); cc++)
1351 {
1352 XdmfSet* xmfSet = xmfGrid->GetSets(cc);
1353 int ghost_value = xmfSet->GetGhost();
1354 if (ghost_value != 0)
1355 {
1356 // skip ghost-sets.
1357 continue;
1358 }
1359 number_of_sets++;
1360 }
1361 if (number_of_sets == 0)
1362 {
1363 return nullptr;
1364 }
1365
1366 vtkMultiBlockDataSet* mb = vtkMultiBlockDataSet::New();
1367 mb->SetNumberOfBlocks(1 + number_of_sets);
1368 mb->SetBlock(0, dataSet);
1369 mb->GetMetaData(static_cast<unsigned int>(0))->Set(vtkCompositeDataSet::NAME(), "Data");
1370
1371 unsigned int current_set_index = 1;
1372 for (int cc = 0; cc < xmfGrid->GetNumberOfSets(); cc++)
1373 {
1374 XdmfSet* xmfSet = xmfGrid->GetSets(cc);
1375 int ghost_value = xmfSet->GetGhost();
1376 if (ghost_value != 0)
1377 {
1378 // skip ghost-sets.
1379 continue;
1380 }
1381
1382 const char* setName = xmfSet->GetName();
1383 mb->GetMetaData(current_set_index)->Set(vtkCompositeDataSet::NAME(), setName);
1384 if (!this->Domain->GetSetsSelection()->ArrayIsEnabled(setName))
1385 {
1386 continue;
1387 }
1388
1389 // Okay now we have an enabled set. Create a new dataset for it
1390 vtkDataSet* set = nullptr;
1391
1392 XdmfInt32 setType = xmfSet->GetSetType();
1393 switch (setType)
1394 {
1395 case XDMF_SET_TYPE_NODE:
1396 set = this->ExtractPoints(xmfSet, dataSet);
1397 break;
1398
1399 case XDMF_SET_TYPE_CELL:
1400 set = this->ExtractCells(xmfSet, dataSet);
1401 break;
1402
1403 case XDMF_SET_TYPE_FACE:
1404 set = this->ExtractFaces(xmfSet, dataSet);
1405 break;
1406
1407 case XDMF_SET_TYPE_EDGE:
1408 set = this->ExtractEdges(xmfSet, dataSet);
1409 break;
1410 }
1411
1412 if (set)
1413 {
1414 mb->SetBlock(current_set_index, set);
1415 set->Delete();
1416 }
1417 current_set_index++;
1418 }
1419 return mb;
1420 }
1421
1422 //------------------------------------------------------------------------------
ExtractPoints(XdmfSet * xmfSet,vtkDataSet * dataSet)1423 vtkDataSet* vtkXdmfHeavyData::ExtractPoints(XdmfSet* xmfSet, vtkDataSet* dataSet)
1424 {
1425 // TODO: How to handle structured datasets with update_extents/strides etc.
1426 // Do they too always produce vtkUniformGrid or do we want to produce
1427 // structured dataset
1428
1429 // Read heavy data. We cannot do anything smart if update_extents or stride
1430 // is specified here. We have to read the entire set and then prune it.
1431 xmfSet->Update();
1432
1433 XdmfArray* xmfIds = xmfSet->GetIds();
1434 XdmfInt64 numIds = xmfIds->GetNumberOfElements();
1435 std::vector<XdmfInt64> ids(numIds + 1);
1436 xmfIds->GetValues(0, ids.data(), numIds);
1437
1438 // release heavy data.
1439 xmfSet->Release();
1440
1441 vtkUnstructuredGrid* output = vtkUnstructuredGrid::New();
1442 vtkPoints* outputPoints = vtkPoints::New();
1443 outputPoints->SetNumberOfPoints(numIds);
1444 output->SetPoints(outputPoints);
1445 outputPoints->Delete();
1446
1447 vtkIdType numInPoints = dataSet->GetNumberOfPoints();
1448 for (vtkIdType kk = 0; kk < numIds; kk++)
1449 {
1450 if (ids[kk] < 0 || ids[kk] > numInPoints)
1451 {
1452 vtkWarningWithObjectMacro(this->Reader, "No such cell or point exists: " << ids[kk]);
1453 continue;
1454 }
1455 double point_location[3];
1456 dataSet->GetPoint(ids[kk], point_location);
1457 outputPoints->SetPoint(kk, point_location);
1458 }
1459 ids.clear(); // done with ids
1460
1461 // Read node-centered attributes that may be defined on this set.
1462 int numAttributes = xmfSet->GetNumberOfAttributes();
1463 for (int cc = 0; cc < numAttributes; cc++)
1464 {
1465 XdmfAttribute* xmfAttribute = xmfSet->GetAttribute(cc);
1466 const char* attrName = xmfAttribute->GetName();
1467 int attrCenter = xmfAttribute->GetAttributeCenter();
1468 if (attrCenter != XDMF_ATTRIBUTE_CENTER_NODE)
1469 {
1470 continue;
1471 }
1472 vtkDataArray* array = this->ReadAttribute(xmfAttribute, 1, nullptr);
1473 if (array)
1474 {
1475 array->SetName(attrName);
1476 output->GetPointData()->AddArray(array);
1477 array->Delete();
1478 }
1479 }
1480
1481 std::vector<vtkIdType> vtk_cell_ids(numIds);
1482 std::iota(vtk_cell_ids.begin(), vtk_cell_ids.end(), 0);
1483 output->InsertNextCell(VTK_POLY_VERTEX, numIds, vtk_cell_ids.data());
1484
1485 return output;
1486 }
1487
1488 //------------------------------------------------------------------------------
ExtractCells(XdmfSet * xmfSet,vtkDataSet * dataSet)1489 vtkDataSet* vtkXdmfHeavyData::ExtractCells(XdmfSet* xmfSet, vtkDataSet* dataSet)
1490 {
1491 // TODO: How to handle structured datasets with update_extents/strides etc.
1492 // Do they too always produce vtkUniformGrid or do we want to produce
1493 // structured dataset
1494
1495 // Read heavy data.
1496 xmfSet->Update();
1497
1498 XdmfArray* xmfIds = xmfSet->GetIds();
1499 XdmfInt64 numIds = xmfIds->GetNumberOfElements();
1500
1501 vtkIdTypeArray* ids = vtkIdTypeArray::New();
1502 ids->SetNumberOfComponents(1);
1503 ids->SetNumberOfTuples(numIds);
1504 xmfIds->GetValues(0, (vtkXdmfIdType*)ids->GetPointer(0), numIds);
1505
1506 // release heavy data.
1507 xmfSet->Release();
1508
1509 // We directly use vtkExtractSelectedIds for extract cells since the logic to
1510 // extract cells it no trivial (like extracting points).
1511 vtkSelectionNode* selNode = vtkSelectionNode::New();
1512 selNode->SetContentType(vtkSelectionNode::INDICES);
1513 selNode->SetFieldType(vtkSelectionNode::CELL);
1514 selNode->SetSelectionList(ids);
1515
1516 vtkSelection* sel = vtkSelection::New();
1517 sel->AddNode(selNode);
1518 selNode->Delete();
1519
1520 vtkExtractSelectedIds* extractCells = vtkExtractSelectedIds::New();
1521 extractCells->SetInputData(0, dataSet);
1522 extractCells->SetInputData(1, sel);
1523 extractCells->Update();
1524
1525 vtkDataSet* output = vtkDataSet::SafeDownCast(extractCells->GetOutput()->NewInstance());
1526 output->CopyStructure(vtkDataSet::SafeDownCast(extractCells->GetOutput()));
1527
1528 sel->Delete();
1529 extractCells->Delete();
1530 ids->Delete();
1531
1532 // Read cell-centered attributes that may be defined on this set.
1533 int numAttributes = xmfSet->GetNumberOfAttributes();
1534 for (int cc = 0; cc < numAttributes; cc++)
1535 {
1536 XdmfAttribute* xmfAttribute = xmfSet->GetAttribute(cc);
1537 const char* attrName = xmfAttribute->GetName();
1538 int attrCenter = xmfAttribute->GetAttributeCenter();
1539 if (attrCenter != XDMF_ATTRIBUTE_CENTER_CELL)
1540 {
1541 continue;
1542 }
1543 vtkDataArray* array = this->ReadAttribute(xmfAttribute, 1, nullptr);
1544 if (array)
1545 {
1546 array->SetName(attrName);
1547 output->GetCellData()->AddArray(array);
1548 array->Delete();
1549 }
1550 }
1551
1552 return output;
1553 }
1554
1555 //------------------------------------------------------------------------------
ExtractFaces(XdmfSet * xmfSet,vtkDataSet * dataSet)1556 vtkDataSet* vtkXdmfHeavyData::ExtractFaces(XdmfSet* xmfSet, vtkDataSet* dataSet)
1557 {
1558 xmfSet->Update();
1559
1560 XdmfArray* xmfIds = xmfSet->GetIds();
1561 XdmfArray* xmfCellIds = xmfSet->GetCellIds();
1562
1563 XdmfInt64 numFaces = xmfIds->GetNumberOfElements();
1564
1565 // ids is a 2 component array were each tuple is (cell-id, face-id).
1566 vtkIdTypeArray* ids = vtkIdTypeArray::New();
1567 ids->SetNumberOfComponents(2);
1568 ids->SetNumberOfTuples(numFaces);
1569 xmfCellIds->GetValues(0, (vtkXdmfIdType*)ids->GetPointer(0), numFaces, 1, 2);
1570 xmfIds->GetValues(0, (vtkXdmfIdType*)ids->GetPointer(1), numFaces, 1, 2);
1571
1572 vtkPolyData* output = vtkPolyData::New();
1573 vtkCellArray* polys = vtkCellArray::New();
1574 output->SetPolys(polys);
1575 polys->Delete();
1576
1577 vtkPoints* outPoints = vtkPoints::New();
1578 output->SetPoints(outPoints);
1579 outPoints->Delete();
1580
1581 vtkMergePoints* mergePoints = vtkMergePoints::New();
1582 mergePoints->InitPointInsertion(outPoints, dataSet->GetBounds());
1583
1584 for (vtkIdType cc = 0; cc < numFaces; cc++)
1585 {
1586 vtkIdType cellId = ids->GetValue(cc * 2);
1587 vtkIdType faceId = ids->GetValue(cc * 2 + 1);
1588 vtkCell* cell = dataSet->GetCell(cellId);
1589 if (!cell)
1590 {
1591 vtkWarningWithObjectMacro(this->Reader, "Invalid cellId: " << cellId);
1592 continue;
1593 }
1594 vtkCell* face = cell->GetFace(faceId);
1595 if (!face)
1596 {
1597 vtkWarningWithObjectMacro(this->Reader, "Invalid faceId " << faceId << " on cell " << cellId);
1598 continue;
1599 }
1600
1601 // Now insert this face a new cell in the output dataset.
1602 vtkIdType numPoints = face->GetNumberOfPoints();
1603 vtkPoints* facePoints = face->GetPoints();
1604 std::vector<vtkIdType> outputPts(numPoints + 1);
1605 for (vtkIdType kk = 0; kk < numPoints; kk++)
1606 {
1607 mergePoints->InsertUniquePoint(facePoints->GetPoint(kk), outputPts[kk]);
1608 }
1609 polys->InsertNextCell(numPoints, outputPts.data());
1610 }
1611
1612 ids->Delete();
1613 xmfSet->Release();
1614 mergePoints->Delete();
1615
1616 // Read face-centered attributes that may be defined on this set.
1617 int numAttributes = xmfSet->GetNumberOfAttributes();
1618 for (int cc = 0; cc < numAttributes; cc++)
1619 {
1620 XdmfAttribute* xmfAttribute = xmfSet->GetAttribute(cc);
1621 const char* attrName = xmfAttribute->GetName();
1622 int attrCenter = xmfAttribute->GetAttributeCenter();
1623 if (attrCenter != XDMF_ATTRIBUTE_CENTER_FACE)
1624 {
1625 continue;
1626 }
1627 vtkDataArray* array = this->ReadAttribute(xmfAttribute, 1, nullptr);
1628 if (array)
1629 {
1630 array->SetName(attrName);
1631 output->GetCellData()->AddArray(array);
1632 array->Delete();
1633 }
1634 }
1635
1636 return output;
1637 }
1638
1639 //------------------------------------------------------------------------------
ExtractEdges(XdmfSet * xmfSet,vtkDataSet * dataSet)1640 vtkDataSet* vtkXdmfHeavyData::ExtractEdges(XdmfSet* xmfSet, vtkDataSet* dataSet)
1641 {
1642 xmfSet->Update();
1643
1644 XdmfArray* xmfIds = xmfSet->GetIds();
1645 XdmfArray* xmfCellIds = xmfSet->GetCellIds();
1646 XdmfArray* xmfFaceIds = xmfSet->GetFaceIds();
1647
1648 XdmfInt64 numEdges = xmfIds->GetNumberOfElements();
1649
1650 // ids is a 3 component array were each tuple is (cell-id, face-id, edge-id).
1651 vtkIdTypeArray* ids = vtkIdTypeArray::New();
1652 ids->SetNumberOfComponents(3);
1653 ids->SetNumberOfTuples(numEdges);
1654 xmfCellIds->GetValues(0, (vtkXdmfIdType*)ids->GetPointer(0), numEdges, 1, 3);
1655 xmfFaceIds->GetValues(0, (vtkXdmfIdType*)ids->GetPointer(1), numEdges, 1, 3);
1656 xmfIds->GetValues(0, (vtkXdmfIdType*)ids->GetPointer(2), numEdges, 1, 3);
1657
1658 vtkPolyData* output = vtkPolyData::New();
1659 vtkCellArray* lines = vtkCellArray::New();
1660 output->SetLines(lines);
1661 lines->Delete();
1662
1663 vtkPoints* outPoints = vtkPoints::New();
1664 output->SetPoints(outPoints);
1665 outPoints->Delete();
1666
1667 vtkMergePoints* mergePoints = vtkMergePoints::New();
1668 mergePoints->InitPointInsertion(outPoints, dataSet->GetBounds());
1669
1670 for (vtkIdType cc = 0; cc < numEdges; cc++)
1671 {
1672 vtkIdType cellId = ids->GetValue(cc * 3);
1673 vtkIdType faceId = ids->GetValue(cc * 3 + 1);
1674 vtkIdType edgeId = ids->GetValue(cc * 3 + 2);
1675 vtkCell* cell = dataSet->GetCell(cellId);
1676 if (!cell)
1677 {
1678 vtkWarningWithObjectMacro(this->Reader, "Invalid cellId: " << cellId);
1679 continue;
1680 }
1681 vtkCell* face = cell->GetFace(faceId);
1682 if (!face)
1683 {
1684 vtkWarningWithObjectMacro(this->Reader, "Invalid faceId " << faceId << " on cell " << cellId);
1685 continue;
1686 }
1687 vtkCell* edge = cell->GetEdge(edgeId);
1688 if (!edge)
1689 {
1690 vtkWarningWithObjectMacro(this->Reader,
1691 "Invalid edgeId " << edgeId << " on face " << faceId << " on cell " << cellId);
1692 continue;
1693 }
1694
1695 // Now insert this edge as a new cell in the output dataset.
1696 vtkIdType numPoints = edge->GetNumberOfPoints();
1697 vtkPoints* edgePoints = edge->GetPoints();
1698 std::vector<vtkIdType> outputPts(numPoints + 1);
1699 for (vtkIdType kk = 0; kk < numPoints; kk++)
1700 {
1701 mergePoints->InsertUniquePoint(edgePoints->GetPoint(kk), outputPts[kk]);
1702 }
1703 lines->InsertNextCell(numPoints, outputPts.data());
1704 }
1705
1706 ids->Delete();
1707 xmfSet->Release();
1708 mergePoints->Delete();
1709
1710 // Read edge-centered attributes that may be defined on this set.
1711 int numAttributes = xmfSet->GetNumberOfAttributes();
1712 for (int cc = 0; cc < numAttributes; cc++)
1713 {
1714 XdmfAttribute* xmfAttribute = xmfSet->GetAttribute(cc);
1715 const char* attrName = xmfAttribute->GetName();
1716 int attrCenter = xmfAttribute->GetAttributeCenter();
1717 if (attrCenter != XDMF_ATTRIBUTE_CENTER_EDGE)
1718 {
1719 continue;
1720 }
1721 vtkDataArray* array = this->ReadAttribute(xmfAttribute, 1, nullptr);
1722 if (array)
1723 {
1724 array->SetName(attrName);
1725 output->GetCellData()->AddArray(array);
1726 array->Delete();
1727 }
1728 }
1729
1730 return output;
1731 }
1732