// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen // SPDX-License-Identifier: BSD-3-Clause #include "vtkHDFReaderImplementation.h" #include #include #include "vtkAMRBox.h" #include "vtkAMRUtilities.h" #include "vtkBitArray.h" #include "vtkCellData.h" #include "vtkDataArrayRange.h" #include "vtkDataArraySelection.h" #include "vtkDataAssembly.h" #include "vtkDataObject.h" #include "vtkFieldData.h" #include "vtkHDF5ScopedHandle.h" #include "vtkHDFUtilities.h" #include "vtkHyperTree.h" #include "vtkHyperTreeGrid.h" #include "vtkHyperTreeGridNonOrientedCursor.h" #include "vtkImageData.h" #include "vtkMemoryResourceStream.h" #include "vtkMultiBlockDataSet.h" #include "vtkOverlappingAMR.h" #include "vtkPartitionedDataSet.h" #include "vtkPartitionedDataSetCollection.h" #include "vtkPolyData.h" #include "vtkStringFormatter.h" #include "vtkUniformGrid.h" #include "vtkUnstructuredGrid.h" #include #include #include #include //------------------------------------------------------------------------------ VTK_ABI_NAMESPACE_BEGIN namespace { constexpr vtkIdType BYTE_SIZE = 8; auto bitSizeToBytes = [](vtkIdType numBits) { return (numBits + BYTE_SIZE - 1) / BYTE_SIZE; // Integer 'ceil' }; } //------------------------------------------------------------------------------ vtkHDFReader::Implementation::Implementation(vtkHDFReader* reader) : File(-1) , VTKGroup(-1) , DataSetType(-1) , NumberOfPieces(-1) , Reader(reader) { std::fill(this->AttributeDataGroup.begin(), this->AttributeDataGroup.end(), -1); std::fill(this->Version.begin(), this->Version.end(), 0); } //------------------------------------------------------------------------------ vtkHDFReader::Implementation::~Implementation() { this->Close(); } //------------------------------------------------------------------------------ std::vector vtkHDFReader::Implementation::GetDimensions(const char* datasetName) { return vtkHDFUtilities::GetDimensions(this->File, datasetName); } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::Open(const char* fileName) { if (!fileName) { vtkErrorWithObjectMacro(this->Reader, "Invalid filename: " << fileName); return false; } if (this->FileName.empty() || this->FileName != fileName || this->File < 0) { this->FileName = fileName; if (this->File >= 0) { this->Close(); } if (!vtkHDFUtilities::Open(fileName, this->File)) { return false; } return this->RetrieveHDFInformation(vtkHDFUtilities::VTKHDF_ROOT_PATH); } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::Open(vtkResourceStream* stream) { if (!stream) { vtkErrorWithObjectMacro(this->Reader, "Stream is nullptr"); return false; } if (!this->Stream || this->Stream != stream || this->File < 0) { this->Stream = stream; if (this->File >= 0) { this->Close(); } auto memStream = vtkMemoryResourceStream::SafeDownCast(stream); if (!memStream) { // Copy to a mem stream if needed stream->Seek(0, vtkResourceStream::SeekDirection::End); std::size_t size = stream->Tell(); stream->Seek(0, vtkResourceStream::SeekDirection::Begin); std::vector tempBuffer; tempBuffer.resize(size); stream->Read(tempBuffer.data(), size); this->LocalMemStream->SetBuffer(std::move(tempBuffer)); memStream = this->LocalMemStream; } if (!vtkHDFUtilities::Open(memStream, this->File)) { return false; } return this->RetrieveHDFInformation(vtkHDFUtilities::VTKHDF_ROOT_PATH); } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::OpenGroupAsVTKGroup(const std::string& groupPath) { if (this->VTKGroup >= 0) { H5Gclose(this->VTKGroup); this->VTKGroup = -1; } if ((this->VTKGroup = H5Gopen(this->File, groupPath.c_str(), H5P_DEFAULT)) < 0) { // the file doesn't exist or we try to read a non-VTKHDF file return false; } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::RetrieveHDFInformation(const std::string& rootName) { this->CloseMemberGroups(); return vtkHDFUtilities::RetrieveHDFInformation(this->File, this->VTKGroup, rootName, this->Version, this->DataSetType, this->NumberOfPieces, this->AttributeDataGroup); } //------------------------------------------------------------------------------ std::size_t vtkHDFReader::Implementation::GetNumberOfSteps() { if (this->File < 0) { vtkErrorWithObjectMacro(this->Reader, "Cannot get number of steps if the file is not open"); return 0; } return vtkHDFUtilities::GetNumberOfSteps(this->VTKGroup); } //------------------------------------------------------------------------------ int vtkHDFReader::Implementation::GetNumberOfPieces(vtkIdType step) { if (step < 0 || this->GetNumberOfSteps() == 1 || H5Lexists(this->VTKGroup, "Steps/NumberOfParts", H5P_DEFAULT) <= 0) { return this->NumberOfPieces; } std::vector buffer = this->GetMetadata("Steps/NumberOfParts", 1, step); if (buffer.empty()) { vtkErrorWithObjectMacro( nullptr, "Could not read step " << step << " in NumberOfParts data set."); return -1; } this->NumberOfPieces = buffer[0]; return this->NumberOfPieces; } //------------------------------------------------------------------------------ void vtkHDFReader::Implementation::Close() { this->DataSetType = -1; this->NumberOfPieces = 0; this->CloseMemberGroups(); std::fill(this->Version.begin(), this->Version.end(), 0); if (this->File >= 0) { H5Fclose(this->File); this->File = -1; } } //------------------------------------------------------------------------------ void vtkHDFReader::Implementation::CloseMemberGroups() { if (this->VTKGroup >= 0) { H5Gclose(this->VTKGroup); this->VTKGroup = -1; } for (size_t i = 0; i < this->AttributeDataGroup.size(); ++i) { if (this->AttributeDataGroup[i] >= 0) { H5Gclose(this->AttributeDataGroup[i]); this->AttributeDataGroup[i] = -1; } } } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::GetPartitionExtent(hsize_t partitionIndex, int* extent) { constexpr int RANK = 2; const char* datasetName = "/VTKHDF/Extents"; // create the memory space hsize_t dimsm[RANK]; dimsm[0] = 1; dimsm[1] = 6; vtkHDF::ScopedH5SHandle memspace = H5Screate_simple(RANK, dimsm, nullptr); if (memspace < 0) { vtkErrorWithObjectMacro(this->Reader, << "Error H5Screate_simple for memory space"); return false; } // create the file dataspace + hyperslab vtkHDF::ScopedH5DHandle dataset = H5Dopen(this->File, datasetName, H5P_DEFAULT); if (dataset < 0) { vtkErrorWithObjectMacro(this->Reader, << std::string("Cannot open ") + datasetName); return false; } hsize_t start[RANK] = { partitionIndex, 0 }, count[RANK] = { 1, 2 }; vtkHDF::ScopedH5SHandle dataspace = H5Dget_space(dataset); if (dataspace < 0) { vtkErrorWithObjectMacro( this->Reader, << std::string("Cannot get space for dataset ") + datasetName); return false; } if (H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, start, nullptr, count, nullptr) < 0) { vtkErrorWithObjectMacro( this->Reader, << std::string("Error selecting hyperslab for ") + datasetName); return false; } // read hyperslab if (H5Dread(dataset, H5T_NATIVE_INT, memspace, dataspace, H5P_DEFAULT, extent) < 0) { vtkErrorWithObjectMacro( this->Reader, << std::string("Error reading hyperslab from ") + datasetName); return false; } return true; } //------------------------------------------------------------------------------ template bool vtkHDFReader::Implementation::GetAttribute( const char* attributeName, size_t numberOfElements, T* value) { return vtkHDFUtilities::GetAttribute(this->VTKGroup, attributeName, numberOfElements, value); } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::HasAttribute(const char* groupName, const char* attributeName) { vtkHDF::ScopedH5GHandle groupID = H5Gopen(this->File, groupName, H5P_DEFAULT); if (groupID < 0) { return false; } return H5Aexists(groupID, attributeName) > 0; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::HasDataset(const char* datasetName) { return H5Lexists(this->VTKGroup, datasetName, H5P_DEFAULT) > 0; } //------------------------------------------------------------------------------ std::vector vtkHDFReader::Implementation::GetArrayNames(int attributeType) { return vtkHDFUtilities::GetArrayNames(this->AttributeDataGroup, attributeType); } //------------------------------------------------------------------------------ std::vector vtkHDFReader::Implementation::GetOrderedChildrenOfGroup( const std::string& path) { return vtkHDFUtilities::GetOrderedChildrenOfGroup(this->VTKGroup, path); } //------------------------------------------------------------------------------ vtkDataArray* vtkHDFReader::Implementation::NewArray( int attributeType, const char* name, const std::vector& fileExtent) { return vtkHDFUtilities::NewArrayForGroup( this->AttributeDataGroup[attributeType], name, fileExtent); } //------------------------------------------------------------------------------ vtkDataArray* vtkHDFReader::Implementation::NewArray( int attributeType, const char* name, hsize_t offset, hsize_t size) { std::vector fileExtent = { offset, offset + size }; return vtkHDFUtilities::NewArrayForGroup( this->AttributeDataGroup[attributeType], name, fileExtent); } //------------------------------------------------------------------------------ vtkAbstractArray* vtkHDFReader::Implementation::NewFieldArray( const char* name, vtkIdType offset, vtkIdType size, vtkIdType dimMaxSize) { return vtkHDFUtilities::NewFieldArray(this->AttributeDataGroup, name, offset, size, dimMaxSize); } //------------------------------------------------------------------------------ vtkDataArray* vtkHDFReader::Implementation::NewMetadataArray( const char* name, hsize_t offset, hsize_t size) { std::vector fileExtent = { offset, offset + size }; return vtkHDFUtilities::NewArrayForGroup(this->VTKGroup, name, fileExtent); } //------------------------------------------------------------------------------ std::vector vtkHDFReader::Implementation::GetMetadata( const char* name, hsize_t size, hsize_t offset) { return vtkHDFUtilities::GetMetadata(this->VTKGroup, name, size, offset); } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::IsPathSoftLink(const std::string& path) { H5L_info_t object; auto err = H5Lget_info(this->File, path.c_str(), &object, H5P_DEFAULT); if (err < 0) { vtkWarningWithObjectMacro(this->Reader, "Can't open '" << path << "' link."); return false; } return object.type == H5L_TYPE_SOFT; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::FillAssembly(vtkDataAssembly* assembly) { if (this->DataSetType != VTK_PARTITIONED_DATA_SET_COLLECTION && this->DataSetType != VTK_MULTIBLOCK_DATA_SET) { vtkErrorWithObjectMacro(this->Reader, "Wrong data set type, expected " << VTK_PARTITIONED_DATA_SET_COLLECTION << ", but got: " << this->DataSetType); return false; } std::string assemblyPath = vtkHDFUtilities::VTKHDF_ROOT_PATH + "/Assembly"; vtkHDF::ScopedH5GHandle assemblyID = H5Gopen(this->VTKGroup, assemblyPath.c_str(), H5P_DEFAULT); if (assemblyID <= H5I_INVALID_HID) { vtkErrorWithObjectMacro( this->Reader, "Can't open 'Assembly' group. A valid Composite vtkHDF file should have it."); return false; } return this->FillAssembly(assembly, this->VTKGroup, 0, assemblyPath); } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::FillAssembly( vtkDataAssembly* assembly, hid_t assemblyHandle, int assemblyID, std::string path) { vtkHDF::ScopedH5GHandle currentHandle = H5Gopen(assemblyHandle, path.c_str(), H5P_DEFAULT); if (currentHandle < H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open '" << path << "' group. A valid Composite vtkHDF file should have it."); return false; } std::vector childrenPath; H5Literate_by_name(currentHandle, path.c_str(), H5_INDEX_CRT_ORDER, H5_ITER_INC, nullptr, vtkHDFUtilities::FileInfoCallBack, &childrenPath, H5P_DEFAULT); if (childrenPath.empty()) { return true; } for (auto& childPath : childrenPath) { std::string childFullPath = path + "/" + childPath; vtkHDF::ScopedH5GHandle childHandle = H5Gopen(currentHandle, childFullPath.c_str(), H5P_DEFAULT); if (childHandle <= H5I_INVALID_HID) { continue; } // Prevent to iterate recursively on the dataset itself by checking if it's a link H5L_info_t object; auto err = H5Lget_info(currentHandle, childFullPath.c_str(), &object, H5P_DEFAULT); if (err < 0) { vtkErrorWithObjectMacro(this->Reader, "Can't open '" << childFullPath << "' link."); return false; } if (object.type == H5L_TYPE_SOFT) { int index = 0; vtkHDFUtilities::GetAttribute(childHandle, "Index", 1, &index); assembly->AddDataSetIndex(assemblyID, index); } else { int groupIndex = assembly->AddNode(childPath.c_str(), assemblyID); this->FillAssembly(assembly, currentHandle, groupIndex, childFullPath); } } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ComputeAMRBlocksPerLevels(unsigned int maxLevel) { this->AMRInformation.Clear(); unsigned int level = 0; while (level < maxLevel) { std::stringstream levelGroupName; levelGroupName << "Level" << level; if (H5Lexists(this->VTKGroup, levelGroupName.str().c_str(), H5P_DEFAULT) <= 0) { // The level does not exist, just exit. return true; } vtkHDF::ScopedH5GHandle levelGroupID = H5Gopen(this->VTKGroup, levelGroupName.str().c_str(), H5P_DEFAULT); if (levelGroupID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open group Level" << level); return false; } if (H5Lexists(levelGroupID, "AMRBox", H5P_DEFAULT) <= 0) { vtkErrorWithObjectMacro(this->Reader, "No AMRBox dataset at Level" << level); return false; } vtkHDF::ScopedH5DHandle amrBoxDatasetID = H5Dopen(levelGroupID, "AMRBox", H5P_DEFAULT); if (amrBoxDatasetID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't find AMRBox dataset at Level" << level); return false; } vtkHDF::ScopedH5SHandle spaceID = H5Dget_space(amrBoxDatasetID); if (spaceID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't get space of AMRBox dataset at Level" << level); return false; } std::array dimensions; if (H5Sget_simple_extent_dims(spaceID, dimensions.data(), nullptr) <= 0) { vtkErrorWithObjectMacro( this->Reader, "Can't get space dimensions of AMRBox dataset at Level" << level); return false; } this->AMRInformation.BlocksPerLevel.push_back(dimensions[0]); ++level; } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::GetImageAttributes( int WholeExtent[6], double Origin[3], double Spacing[3]) { if (!this->GetAttribute("WholeExtent", 6, WholeExtent)) { vtkErrorWithObjectMacro(this->Reader, "Could not get WholeExtent attribute"); return false; } if (!this->GetAttribute("Origin", 3, Origin)) { vtkErrorWithObjectMacro(this->Reader, "Could not get Origin attribute"); return false; } if (!this->GetAttribute("Spacing", 3, Spacing)) { vtkErrorWithObjectMacro(this->Reader, "Could not get Spacing attribute"); return false; } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ComputeAMROffsetsPerLevels( vtkDataArraySelection* dataArraySelection[3], vtkIdType step, unsigned int maxLevel) { if (!dataArraySelection) { return false; } this->AMRInformation.Clear(); if (this->DataSetType != VTK_OVERLAPPING_AMR) { vtkWarningWithObjectMacro( this->Reader, "Bad usage of this method. Should only be use for OverlappingAMR"); return true; } vtkIdType numberOfSteps = this->GetNumberOfSteps(); unsigned int level = 0; while (level < maxLevel) { std::stringstream levelGroupName; levelGroupName << "Steps/Level" << level; if (H5Lexists(this->VTKGroup, levelGroupName.str().c_str(), H5P_DEFAULT) <= 0) { // The level does not exist, just exit. return true; } vtkHDF::ScopedH5GHandle levelGroupID = H5Gopen(this->VTKGroup, levelGroupName.str().c_str(), H5P_DEFAULT); if (levelGroupID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open group Level" << level); return false; } std::vector numberOfBox; numberOfBox.resize(numberOfSteps); std::vector boxOffsets; boxOffsets.resize(numberOfSteps); // Due to an error in the VTKHDF File Format 2.2, few array names miss an 's' at // the end. As we should support any VTKHDF version 2.X, we should handle it. // It concerns: NumberOfAMRBoxes, AMRBoxOffsets and Point/Cell/FieldDataOffsets std::string numberOfBoxesDatasetName = "NumberOfAMRBoxes"; std::string amrboxOffsetsDatasetName = "AMRBoxOffsets"; std::array groupNames = { "PointDataOffsets", "CellDataOffsets", "FieldDataOffsets" }; if (this->Version[0] == 2 && this->Version[1] == 2) { numberOfBoxesDatasetName = "NumberOfAMRBox"; amrboxOffsetsDatasetName = "AMRBoxOffset"; groupNames = { "PointDataOffset", "CellDataOffset", "FieldDataOffset" }; } if (H5Lexists(levelGroupID, numberOfBoxesDatasetName.c_str(), H5P_DEFAULT) <= 0) { vtkErrorWithObjectMacro( this->Reader, "No " << numberOfBoxesDatasetName << " dataset at Steps/Level" << level); return false; } vtkHDF::ScopedH5DHandle amrBoxDatasetID = H5Dopen(levelGroupID, numberOfBoxesDatasetName.c_str(), H5P_DEFAULT); if (amrBoxDatasetID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't find " << numberOfBoxesDatasetName << " dataset at Steps/Level" << level); return false; } if (H5Dread( amrBoxDatasetID, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, numberOfBox.data()) < 0) { vtkErrorWithObjectMacro(this->Reader, << "Error reading hyperslab"); return false; } if (H5Lexists(levelGroupID, amrboxOffsetsDatasetName.c_str(), H5P_DEFAULT) <= 0) { vtkErrorWithObjectMacro( this->Reader, "No " << amrboxOffsetsDatasetName << " dataset at Steps/Level" << level); return false; } vtkHDF::ScopedH5DHandle boxOffsetID = H5Dopen(levelGroupID, amrboxOffsetsDatasetName.c_str(), H5P_DEFAULT); if (boxOffsetID == H5I_INVALID_HID) { vtkErrorWithObjectMacro( this->Reader, "Can't " << amrboxOffsetsDatasetName << " dataset at Steps/Level" << level); return false; } if (H5Dread(boxOffsetID, VTK_ID_H5T, H5S_ALL, H5S_ALL, H5P_DEFAULT, boxOffsets.data()) < 0) { vtkErrorWithObjectMacro(this->Reader, << "Error reading hyperslab from "); return false; } this->AMRInformation.BlocksPerLevel.push_back(numberOfBox[step]); this->AMRInformation.BlockOffsetsPerLevel.push_back((boxOffsets[step])); for (int attributeType = vtkDataObject::AttributeTypes::POINT; attributeType <= vtkDataObject::AttributeTypes::FIELD; ++attributeType) { if (H5Lexists(levelGroupID, groupNames[attributeType], H5P_DEFAULT) <= 0) { // These groups are optional continue; } vtkHDF::ScopedH5GHandle cellOffsetID = H5Gopen(levelGroupID, groupNames[attributeType], H5P_DEFAULT); if (cellOffsetID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open " << groupNames[attributeType] << " group at Steps/Level" << level); return false; } const std::vector arrayNames = this->GetArrayNames(attributeType); for (const std::string& name : arrayNames) { if (!dataArraySelection[attributeType]->ArrayIsEnabled(name.c_str())) { continue; } if (H5Lexists(cellOffsetID, name.c_str(), H5P_DEFAULT) <= 0) { vtkErrorWithObjectMacro( this->Reader, "Can't open " << name << " group at Steps/Level" << level); return false; } vtkHDF::ScopedH5DHandle constantID = H5Dopen(cellOffsetID, name.c_str(), H5P_DEFAULT); if (constantID == H5I_INVALID_HID) { vtkErrorWithObjectMacro( this->Reader, "Can't open " << name << " dataset at Steps/Level" << level); return false; } std::vector cellOffsets; cellOffsets.resize(numberOfSteps); if (H5Dread(constantID, VTK_ID_H5T, H5S_ALL, H5S_ALL, H5P_DEFAULT, cellOffsets.data()) < 0) { vtkErrorWithObjectMacro(this->Reader, << "Error reading hyperslab from "); return false; } switch (attributeType) { case vtkDataObject::AttributeTypes::POINT: this->AMRInformation.PointOffsetsPerLevel[name].push_back(cellOffsets[step]); break; case vtkDataObject::AttributeTypes::CELL: this->AMRInformation.CellOffsetsPerLevel[name].push_back(cellOffsets[step]); break; case vtkDataObject::AttributeTypes::FIELD: this->AMRInformation.FieldOffsetsPerLevel[name].push_back(cellOffsets[step]); if (step + 1 < static_cast(cellOffsets.size())) { vtkIdType fieldSize = cellOffsets[step + 1] - cellOffsets[step]; this->AMRInformation.FieldSizesPerLevel[name].push_back(fieldSize); } else { this->AMRInformation.FieldSizesPerLevel[name].push_back(-1); } break; } } } ++level; } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadLevelTopology(unsigned int level, const std::string& levelGroupName, vtkOverlappingAMR* data, double origin[3], bool isTemporalData) { vtkHDF::ScopedH5GHandle levelGroupID = H5Gopen(this->VTKGroup, levelGroupName.c_str(), H5P_DEFAULT); if (levelGroupID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open group Level" << level); return false; } std::array spacing = { 0, 0, 0 }; if (!this->ReadLevelSpacing(levelGroupID, spacing.data())) { vtkErrorWithObjectMacro( this->Reader, "Error while reading spacing attribute at level " << level); return false; } data->SetSpacing(level, spacing.data()); std::vector amrBoxRawData; if (!this->ReadAMRBoxRawValues(levelGroupID, amrBoxRawData, level, isTemporalData)) { vtkErrorWithObjectMacro(this->Reader, "Error while reading AMRBox at level " << level); return false; } if (amrBoxRawData.size() % 6 != 0) { vtkErrorWithObjectMacro(this->Reader, "The size of the \"AMRBox\" dataset at Level" << level << " is not a multiple of 6."); return false; } vtkIdType numberOfDatasets = static_cast(amrBoxRawData.size() / 6); for (vtkIdType dataSetIndex = 0; dataSetIndex < numberOfDatasets; ++dataSetIndex) { int* currentAMRBoxRawData = amrBoxRawData.data() + (6 * dataSetIndex); vtkAMRBox amrBox(currentAMRBoxRawData); data->SetAMRBox(level, dataSetIndex, amrBox); vtkNew dataSet; dataSet->Initialize(); const int* lowCorner = amrBox.GetLoCorner(); double dataSetOrigin[3] = { origin[0] + lowCorner[0] * spacing[0], origin[1] + lowCorner[1] * spacing[1], origin[2] + lowCorner[2] * spacing[2] }; dataSet->SetOrigin(dataSetOrigin); dataSet->SetSpacing(spacing.data()); std::array numberOfNodes = { 0, 0, 0 }; amrBox.GetNumberOfNodes(numberOfNodes.data()); dataSet->SetDimensions(numberOfNodes.data()); data->SetDataSet(level, dataSetIndex, dataSet); } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadLevelData(unsigned int level, const std::string& levelGroupName, vtkOverlappingAMR* data, vtkDataArraySelection* dataArraySelection[3], bool isTemporalData) { vtkHDF::ScopedH5GHandle levelGroupID = H5Gopen(this->VTKGroup, levelGroupName.c_str(), H5P_DEFAULT); if (levelGroupID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open group Level" << level); return false; } // Now read actual data - one array at a time std::array groupNames = { "PointData", "CellData", "FieldData" }; for (int attributeType = vtkDataObject::AttributeTypes::POINT; attributeType <= vtkDataObject::AttributeTypes::FIELD; ++attributeType) { vtkHDF::ScopedH5GHandle groupID = H5Gopen(levelGroupID, groupNames[attributeType], H5P_DEFAULT); if (groupID == H5I_INVALID_HID) { // It's OK to not have groups in the file if there are no data arrays // for that attribute type. continue; } const std::vector arrayNames = this->GetArrayNames(attributeType); for (const std::string& name : arrayNames) { if (!dataArraySelection[attributeType]->ArrayIsEnabled(name.c_str())) { continue; } // Open dataset hid_t tempNativeType = H5I_INVALID_HID; std::vector dims; vtkHDF::ScopedH5DHandle datasetID = vtkHDFUtilities::OpenDataSet(groupID, name.c_str(), &tempNativeType, dims); vtkHDF::ScopedH5THandle nativeType = tempNativeType; if (datasetID < 0) { continue; } // Iterate over all datasets, read data and assign attribute hsize_t dataOffset = 0; hsize_t dataSize = 0; unsigned int numberOfDatasets = data->GetNumberOfBlocks(level); for (unsigned int dataSetIndex = 0; dataSetIndex < numberOfDatasets; ++dataSetIndex) { const vtkAMRBox& amrBox = data->GetAMRBox(level, dataSetIndex); vtkImageData* dataSet = data->GetDataSetAsImageData(level, dataSetIndex); if (dataSet == nullptr) { vtkErrorWithObjectMacro(this->Reader, "Error fetching dataset at level " << level << " and dataSetIndex " << dataSetIndex); return false; } // Here dataSize is the size of the previous dataset read. Offset // is incremented, and a new size is specified after the increment. // This allows for reading AMR's where the size of the blocks vary // inside each level. dataOffset += dataSize; hsize_t cellOffset = 0; switch (attributeType) { case vtkDataObject::AttributeTypes::POINT: dataSize = amrBox.GetNumberOfNodes(); if (isTemporalData) { if (this->AMRInformation.PointOffsetsPerLevel.count(name) == 1) { cellOffset = this->AMRInformation.PointOffsetsPerLevel[name][level]; } } break; case vtkDataObject::AttributeTypes::CELL: dataSize = amrBox.GetNumberOfCells(); if (isTemporalData) { if (this->AMRInformation.CellOffsetsPerLevel.count(name) == 1) { cellOffset = this->AMRInformation.CellOffsetsPerLevel[name][level]; } } break; case vtkDataObject::AttributeTypes::FIELD: dataSize = dims[0] / numberOfDatasets; if (isTemporalData) { if (this->AMRInformation.FieldOffsetsPerLevel.count(name) == 1) { cellOffset = this->AMRInformation.FieldOffsetsPerLevel[name][level] / numberOfDatasets; vtkIdType fieldSize = this->AMRInformation.FieldSizesPerLevel[name][level]; if (fieldSize == -1) { dataSize = (dataSize - cellOffset); } else { dataSize = fieldSize; } dataSize /= numberOfDatasets; } } } std::vector fileExtent = { cellOffset + dataOffset, cellOffset + dataOffset + dataSize }; vtkSmartPointer array; if ((array = vtk::TakeSmartPointer(vtkHDFUtilities::NewArrayForGroup( datasetID, nativeType, dims, fileExtent))) == nullptr) { vtkErrorWithObjectMacro(this->Reader, "Error reading array " << name); return false; } array->SetName(name.c_str()); dataSet->GetAttributesAsFieldData(attributeType)->AddArray(array); } } } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadLevelSpacing(hid_t levelGroupID, double* spacing) { if (!H5Aexists(levelGroupID, "Spacing")) { vtkErrorWithObjectMacro(this->Reader, "\"Spacing\" attribute does not exist."); return false; } vtkHDF::ScopedH5AHandle spacingAttributeID = H5Aopen_name(levelGroupID, "Spacing"); if (spacingAttributeID < 0) { vtkErrorWithObjectMacro(this->Reader, "Can't open \"Spacing\" attribute."); return false; } if (H5Aread(spacingAttributeID, H5T_NATIVE_DOUBLE, spacing) < 0) { vtkErrorWithObjectMacro(this->Reader, "Can't read \"Spacing\" attribute."); return false; } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadAMRBoxRawValues( hid_t levelGroupID, std::vector& amrBoxRawData, int level, bool isTemporalData) { hsize_t startBlock = 0; if (isTemporalData) { startBlock = static_cast(this->AMRInformation.BlockOffsetsPerLevel[level]); } hsize_t numberOfBlock = static_cast(this->AMRInformation.BlocksPerLevel[level]); if (H5Lexists(levelGroupID, "AMRBox", H5P_DEFAULT) <= 0) { vtkErrorWithObjectMacro(this->Reader, "No AMRBox dataset"); return false; } vtkHDF::ScopedH5DHandle amrBoxDatasetID = H5Dopen(levelGroupID, "AMRBox", H5P_DEFAULT); if (amrBoxDatasetID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't open AMRBox dataset"); return false; } vtkHDF::ScopedH5SHandle spaceID = H5Dget_space(amrBoxDatasetID); if (spaceID == H5I_INVALID_HID) { vtkErrorWithObjectMacro(this->Reader, "Can't get space of AMRBox dataset"); return false; } std::array dimensions; if (H5Sget_simple_extent_dims(spaceID, dimensions.data(), nullptr) <= 0) { vtkErrorWithObjectMacro(this->Reader, "Can't get space dimensions of AMRBox dataset"); return false; } if (dimensions[1] != 6) { vtkErrorWithObjectMacro( this->Reader, "Wrong AMRBox dimension, expected 6, got: " << dimensions[1]); return false; } hsize_t startPosition[2] = { startBlock, 0 }; hsize_t count[2] = { numberOfBlock, 6 }; hid_t filSpace = H5Screate_simple(2, dimensions.data(), nullptr); H5Sselect_hyperslab(filSpace, H5S_SELECT_SET, startPosition, nullptr, count, nullptr); startPosition[0] = 0; hid_t memSpace = H5Screate_simple(2, dimensions.data(), nullptr); H5Sselect_hyperslab(memSpace, H5S_SELECT_SET, startPosition, nullptr, count, nullptr); hsize_t numberOfDatasets = numberOfBlock; amrBoxRawData.resize(numberOfDatasets * 6); if (H5Dread( amrBoxDatasetID, H5T_NATIVE_INT, memSpace, filSpace, H5P_DEFAULT, amrBoxRawData.data()) < 0) { vtkErrorWithObjectMacro(this->Reader, "Can't read AMRBox dataset."); return false; } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadAMRTopology(vtkOverlappingAMR* data, unsigned int level, unsigned int maxLevel, double origin[3], bool isTemporalData) { if (this->AMRInformation.BlocksPerLevel.empty()) { return false; } size_t numberOfLoadedLevels = 0; if (maxLevel == 0) { numberOfLoadedLevels = this->AMRInformation.BlocksPerLevel.size(); } else { numberOfLoadedLevels = std::min(this->AMRInformation.BlocksPerLevel.size(), static_cast(maxLevel)); } std::vector blocksPerLevel; blocksPerLevel.reserve(numberOfLoadedLevels); for (size_t i = 0; i < numberOfLoadedLevels; i++) { blocksPerLevel.emplace_back(this->AMRInformation.BlocksPerLevel[i]); } data->Initialize(blocksPerLevel); data->SetOrigin(origin); data->SetGridDescription(vtkStructuredData::VTK_STRUCTURED_XYZ_GRID); while (level < maxLevel) { std::string levelGroupName = "Level" + vtk::to_string(level); if (H5Lexists(this->VTKGroup, levelGroupName.c_str(), H5P_DEFAULT) <= 0) { break; } else { if (!this->ReadLevelTopology(level, levelGroupName, data, origin, isTemporalData)) { vtkErrorWithObjectMacro(this->Reader, << "Can't read group Level" << level); return false; } ++level; } } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadAMRData(vtkOverlappingAMR* data, unsigned int level, unsigned int maxLevel, vtkDataArraySelection* dataArraySelection[3], bool isTemporalData) { while (level < maxLevel) { std::string levelGroupName = "Level" + vtk::to_string(level); if (H5Lexists(this->VTKGroup, levelGroupName.c_str(), H5P_DEFAULT) <= 0) { break; } else { if (!this->ReadLevelData(level, levelGroupName, data, dataArraySelection, isTemporalData)) { vtkErrorWithObjectMacro(this->Reader, << "Can't fill group Level" << level); return false; } ++level; } } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadHyperTreeGridData(vtkHyperTreeGrid* htg, const vtkDataArraySelection* arraySelection, const vtkIdType cellOffset, const vtkIdType treeIdsOffset, const vtkIdType depthOffset, const vtkIdType descriptorOffset, const vtkIdType maskOffset, const vtkIdType partOffset, const vtkIdType verticesPerDepthOffset, const vtkIdType depthLimit, const vtkIdType step) { htg->Initialize(); if (!this->ReadHyperTreeGridMetaInfo(htg)) { return false; } if (!this->ReadHyperTreeGridDimensions(htg)) { return false; } // Read meta-data for the current piece const vtkIdType treeCount = this->GetMetadata("NumberOfTrees", 1, partOffset)[0]; const vtkIdType depthCount = this->GetMetadata("NumberOfDepths", 1, partOffset)[0]; const vtkIdType cellCount = this->GetMetadata("NumberOfCells", 1, partOffset)[0]; const vtkIdType descriptorCount = this->GetMetadata("DescriptorsSize", 1, partOffset)[0]; // Read only the tree ids for the current piece const std::vector treeIds = this->GetMetadata("TreeIds", treeCount, treeIdsOffset); const std::vector depthPerTree = this->GetMetadata("DepthPerTree", treeCount, depthOffset); const std::vector numberOfCellsPerTreeDepth = this->GetMetadata("NumberOfCellsPerTreeDepth", depthCount, verticesPerDepthOffset); // Read Tree Descriptors const std::vector descriptorExtent = { static_cast(descriptorOffset), static_cast(descriptorOffset + ::bitSizeToBytes(descriptorCount)) }; vtkSmartPointer descriptorsByteArray = vtk::TakeSmartPointer(vtkArrayDownCast( vtkHDFUtilities::NewArrayForGroup(this->VTKGroup, "Descriptors", descriptorExtent))); // Build a bit array from the uint8 array, transferring the memory ownership to the bit array vtkNew descriptor; descriptor->SetArray(descriptorsByteArray->GetPointer(0), descriptorCount, 0); descriptorsByteArray->SetArrayFreeFunction(nullptr); const bool hasMask = H5Lexists(this->VTKGroup, "Mask", H5P_DEFAULT) > 0; if (hasMask) { vtkNew mask; mask->Allocate(cellCount); htg->SetMask(mask); } // Create cell arrays std::vector> cellArrays; if (!this->CreateHyperTreeGridCellArrays(htg, cellArrays, arraySelection, cellCount)) { return false; } vtkNew treeCursor; vtkIdType treeDescriptorOffset = 0; vtkIdType inputCellOffset = 0; vtkIdType outputCellOffset = 0; vtkIdType numCellsPerDepthOffset = 0; // Iterate over trees, building them taking depth limiter into consideration. // Also appends masking and cell data for (const vtkIdType& treeId : treeIds) { // Compute both the total number of cells and the number of cells below the depth limit in the // current tree Number of cells in the current tree = sum of the number of cells in all depths vtkIdType lastDepthSize = 0; vtkIdType lastReadableDepthSize = 0; // Taking depth limit into account vtkIdType totalTreeSize = 0; vtkIdType readableTreeSize = 0; // Taking depth limit into account auto numCellsPerTreeDepthIt = numberOfCellsPerTreeDepth.begin() + numCellsPerDepthOffset; for (vtkIdType depth = 0; depth < depthPerTree[treeId]; depth++, numCellsPerTreeDepthIt++) { lastDepthSize = *numCellsPerTreeDepthIt; if (depth < depthLimit) { lastReadableDepthSize = lastDepthSize; readableTreeSize += lastDepthSize; } totalTreeSize += lastDepthSize; } // The descriptor has 1 bit entry for each cell in the tree, except for the last depth const vtkIdType descriptorSize = totalTreeSize - lastDepthSize; const vtkIdType readableDescriptorSize = readableTreeSize - lastReadableDepthSize; htg->InitializeNonOrientedCursor(treeCursor, treeId, true); treeCursor->SetGlobalIndexStart(outputCellOffset); if (!this->AppendCellDataForHyperTree( cellArrays, cellOffset, inputCellOffset, step, readableTreeSize)) { return false; } if (hasMask && !this->AppendMaskForHyperTree(htg, inputCellOffset, maskOffset, readableTreeSize)) { return false; } // Build tree from descriptor treeCursor->GetTree()->BuildFromBreadthFirstOrderDescriptor( descriptor, readableDescriptorSize, treeDescriptorOffset); // Increment input & output offsets treeDescriptorOffset += descriptorSize; inputCellOffset += totalTreeSize; outputCellOffset += readableTreeSize; numCellsPerDepthOffset += depthPerTree[treeId]; } // We pre-allocated `cellCount` items, but with depth limiter, it is possible than less than that // is used. for (auto& array : cellArrays) { array->Squeeze(); } if (hasMask) { htg->GetMask()->Squeeze(); } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadHyperTreeGridMetaInfo(vtkHyperTreeGrid* htg) { // Read general info about this HTG int branchFactor = 0; if (!this->GetAttribute("BranchFactor", 1, &branchFactor)) { vtkErrorWithObjectMacro(this->Reader, << "Cannot find BranchFactor attribute"); return false; } htg->SetBranchFactor(branchFactor); // Read TransposedRootIndexing attribute int transposedRootIndexing = 0; if (H5Aexists(this->VTKGroup, "TransposedRootIndexing") > 0) { this->GetAttribute("TransposedRootIndexing", 1, &transposedRootIndexing); } htg->SetTransposedRootIndexing(transposedRootIndexing != 0); // Read Interface attributes if they exist bool hasInterfaceIntercepts = H5Aexists(this->VTKGroup, "InterfaceInterceptsName") > 0; bool hasInterfaceNormals = H5Aexists(this->VTKGroup, "InterfaceNormalsName") > 0; if (hasInterfaceIntercepts) { std::string interfaceInterceptsName; vtkHDFUtilities::GetStringAttribute( this->VTKGroup, "InterfaceInterceptsName", interfaceInterceptsName); htg->SetInterfaceInterceptsName(interfaceInterceptsName.c_str()); } if (hasInterfaceNormals) { std::string interfaceNormalsName; vtkHDFUtilities::GetStringAttribute( this->VTKGroup, "InterfaceNormalsName", interfaceNormalsName); htg->SetInterfaceNormalsName(interfaceNormalsName.c_str()); } if (hasInterfaceIntercepts && hasInterfaceNormals) { htg->SetHasInterface(true); } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::ReadHyperTreeGridDimensions(vtkHyperTreeGrid* htg) { std::array dimensions; if (!this->GetAttribute("Dimensions", 1, dimensions.data())) { vtkErrorWithObjectMacro(this->Reader, << "Missing dimensions attribute"); return false; } // Read coordinate arrays std::vector coordinates_extent{ 0, static_cast(dimensions[0]) }; auto XCoordinates = vtk::TakeSmartPointer( vtkHDFUtilities::NewArrayForGroup(this->VTKGroup, "XCoordinates", coordinates_extent)); coordinates_extent[1] = static_cast(dimensions[1]); auto YCoordinates = vtk::TakeSmartPointer( vtkHDFUtilities::NewArrayForGroup(this->VTKGroup, "YCoordinates", coordinates_extent)); coordinates_extent[1] = static_cast(dimensions[2]); auto ZCoordinates = vtk::TakeSmartPointer( vtkHDFUtilities::NewArrayForGroup(this->VTKGroup, "ZCoordinates", coordinates_extent)); if (!XCoordinates || !YCoordinates || !ZCoordinates) { vtkErrorWithObjectMacro(this->Reader, << "Missing coordinates array"); return false; } htg->SetDimensions(dimensions.data()); htg->SetXCoordinates(XCoordinates); htg->SetYCoordinates(YCoordinates); htg->SetZCoordinates(ZCoordinates); return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::CreateHyperTreeGridCellArrays(vtkHyperTreeGrid* htg, std::vector>& cellArrays, const vtkDataArraySelection* arraySelection, vtkIdType cellCount) { constexpr int cellType = vtkDataObject::AttributeTypes::CELL; const std::vector arrayNames = this->GetArrayNames(cellType); for (const std::string& arrayName : arrayNames) { if (arraySelection->ArrayIsEnabled(arrayName.c_str())) { // Read a null extent so we get the correct array type vtkSmartPointer array; std::vector nullExtent{ 0, 0 }; if ((array = vtk::TakeSmartPointer( this->NewArray(cellType, arrayName.c_str(), nullExtent))) == nullptr) { vtkErrorWithObjectMacro(nullptr, "Error reading array " << arrayName); return false; } array->SetName(arrayName.c_str()); array->Allocate(cellCount * array->GetNumberOfComponents()); htg->GetCellData()->AddArray(array); cellArrays.emplace_back(array); } } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::AppendCellDataForHyperTree( std::vector>& cellArrays, vtkIdType cellOffset, vtkIdType inputCellOffset, const vtkIdType step, const vtkIdType readableTreeSize) { constexpr int cellType = vtkDataObject::AttributeTypes::CELL; for (auto& array : cellArrays) { std::string arrayName = array->GetName(); vtkIdType startReadOffset = cellOffset + std::max(this->GetArrayOffset(step, cellType, arrayName), 0) + inputCellOffset; std::vector arrayReadExtent{ static_cast(startReadOffset), static_cast(startReadOffset + readableTreeSize) }; vtkSmartPointer addedArray; if ((addedArray = vtk::TakeSmartPointer( this->NewArray(cellType, arrayName.c_str(), arrayReadExtent))) == nullptr) { vtkErrorWithObjectMacro(nullptr, "Error reading array " << arrayName); return false; } // Copy read array at the end of the data array array->InsertTuples(array->GetNumberOfTuples(), readableTreeSize, 0, addedArray); } return true; } //------------------------------------------------------------------------------ bool vtkHDFReader::Implementation::AppendMaskForHyperTree(vtkHyperTreeGrid* htg, vtkIdType inputCellOffset, const vtkIdType maskOffset, const vtkIdType readableTreeSize) { // Read and append mask vtkIdType maskByteStartOffset = maskOffset + (inputCellOffset / BYTE_SIZE); vtkIdType offsetInStartingByte = inputCellOffset % BYTE_SIZE; // Mask for the tree may not start at the start of a byte vtkIdType maskByteSize = ::bitSizeToBytes(readableTreeSize + offsetInStartingByte); std::vector maskExtent = { static_cast(maskByteStartOffset), static_cast(maskByteStartOffset + maskByteSize) }; // Read only the relevant part of the mask, taking depth limiter into account vtkSmartPointer addedMaskByteArray = vtk::TakeSmartPointer(vtkArrayDownCast( vtkHDFUtilities::NewArrayForGroup(this->VTKGroup, "Mask", maskExtent))); vtkNew addedMask; addedMask->SetArray(vtkArrayDownCast(addedMaskByteArray)->GetPointer(0), maskByteSize * BYTE_SIZE, 0); addedMaskByteArray->SetArrayFreeFunction(nullptr); if (htg->GetMask()) { htg->GetMask()->InsertTuples( htg->GetMask()->GetNumberOfTuples(), readableTreeSize, offsetInStartingByte, addedMask); } return true; } //------------------------------------------------------------------------------ vtkDataArray* vtkHDFReader::Implementation::GetStepValues() { if (this->File < 0) { vtkErrorWithObjectMacro(this->Reader, "Cannot get step values if the file is not open"); } return this->GetStepValues(this->VTKGroup); } //------------------------------------------------------------------------------ vtkDataArray* vtkHDFReader::Implementation::GetStepValues(hid_t group) { if (group < 0) { vtkErrorWithObjectMacro(this->Reader, "Cannot get step values from empty group"); } if (H5Lexists(group, "Steps", H5P_DEFAULT) <= 0) { // Steps group does not exist return nullptr; } // Steps group does exist vtkHDF::ScopedH5GHandle steps = H5Gopen(group, "Steps", H5P_DEFAULT); if (steps < 0) { vtkErrorWithObjectMacro(this->Reader, "Could not open steps group"); return nullptr; } std::vector fileExtent; return vtkHDFUtilities::NewArrayForGroup(steps, "Values", fileExtent); } //------------------------------------------------------------------------------ vtkIdType vtkHDFReader::Implementation::GetArrayOffset( vtkIdType step, int attributeType, std::string name) { return vtkHDFUtilities::GetArrayOffset(this->VTKGroup, step, attributeType, name); } //------------------------------------------------------------------------------ std::array vtkHDFReader::Implementation::GetFieldArraySize( vtkIdType step, std::string name) { return vtkHDFUtilities::GetFieldArraySize(this->VTKGroup, step, name); } //------------------------------------------------------------------------------ // explicit template instantiation template bool vtkHDFReader::Implementation::GetAttribute( const char* attributeName, size_t dim, int* value); template bool vtkHDFReader::Implementation::GetAttribute( const char* attributeName, size_t dim, double* value); VTK_ABI_NAMESPACE_END //------------------------------------------------------------------------------ vtkSmartPointer vtkHDFReader::Implementation::GetNewDataSet( int dataSetType, int numPieces) { constexpr std::array partitionedTypes = { VTK_UNSTRUCTURED_GRID, VTK_POLY_DATA, VTK_HYPER_TREE_GRID }; vtkSmartPointer newOutput = nullptr; if (numPieces > 1 && std::find(partitionedTypes.begin(), partitionedTypes.end(), dataSetType) != partitionedTypes.end()) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_IMAGE_DATA) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_UNSTRUCTURED_GRID) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_POLY_DATA) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_HYPER_TREE_GRID) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_OVERLAPPING_AMR) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_PARTITIONED_DATA_SET_COLLECTION) { newOutput = vtkSmartPointer::New(); } else if (dataSetType == VTK_MULTIBLOCK_DATA_SET) { newOutput = vtkSmartPointer::New(); } else { vtkErrorWithObjectMacro(this->Reader, "HDF dataset type unknown: " << dataSetType); return nullptr; } return newOutput; }