diff --git a/docs/changelog/write-uniform-rectilinear.md b/docs/changelog/write-uniform-rectilinear.md new file mode 100644 index 000000000..1d601060d --- /dev/null +++ b/docs/changelog/write-uniform-rectilinear.md @@ -0,0 +1,13 @@ +# Write uniform and rectilinear grids to legacy VTK files + +As a programming convenience, all `vtkm::cont::DataSet` written by +`vtkm::io::VTKDataSetWriter` were written as a structured grid. Although +technically correct, it changed the structure of the data. This meant that +if you wanted to capture data to run elsewhere, it would run as a different +data type. This was particularly frustrating if the data of that structure +was causing problems and you wanted to debug it. + +Now, `VTKDataSetWriter` checks the type of the `CoordinateSystem` to +determine whether the data should be written out as `STRUCTURED_POINTS` +(i.e. a uniform grid), `RECTILINEAR_GRID`, or `STRUCTURED_GRID` +(curvilinear). diff --git a/vtkm/io/VTKDataSetWriter.cxx b/vtkm/io/VTKDataSetWriter.cxx index 6f3eac1a3..43b3aaffa 100644 --- a/vtkm/io/VTKDataSetWriter.cxx +++ b/vtkm/io/VTKDataSetWriter.cxx @@ -34,6 +34,12 @@ namespace { +template +using ArrayHandleRectilinearCoordinates = + vtkm::cont::ArrayHandleCartesianProduct, + vtkm::cont::ArrayHandle, + vtkm::cont::ArrayHandle>; + struct OutputPointsFunctor { private: @@ -135,6 +141,18 @@ private: std::string* Name; }; +template +void WriteDimensions(std::ostream& out, const vtkm::cont::CellSetStructured& cellSet) +{ + auto pointDimensions = cellSet.GetPointDimensions(); + using VTraits = vtkm::VecTraits; + + out << "DIMENSIONS "; + out << VTraits::GetComponent(pointDimensions, 0) << " "; + out << (DIM > 1 ? VTraits::GetComponent(pointDimensions, 1) : 1) << " "; + out << (DIM > 2 ? VTraits::GetComponent(pointDimensions, 2) : 1) << "\n"; +} + void WritePoints(std::ostream& out, const vtkm::cont::DataSet& dataSet) { ///\todo: support other coordinate systems @@ -302,23 +320,89 @@ void WriteDataSetAsUnstructured(std::ostream& out, WriteExplicitCells(out, cellSet); } +template +void WriteDataSetAsStructuredPoints(std::ostream& out, + const vtkm::cont::ArrayHandleUniformPointCoordinates& points, + const vtkm::cont::CellSetStructured& cellSet) +{ + out << "DATASET STRUCTURED_POINTS\n"; + + WriteDimensions(out, cellSet); + + auto portal = points.ReadPortal(); + auto origin = portal.GetOrigin(); + auto spacing = portal.GetSpacing(); + out << "ORIGIN " << origin[0] << " " << origin[1] << " " << origin[2] << "\n"; + out << "SPACING " << spacing[0] << " " << spacing[1] << " " << spacing[2] << "\n"; +} + +template +void WriteDataSetAsRectilinearGrid(std::ostream& out, + const ArrayHandleRectilinearCoordinates& points, + const vtkm::cont::CellSetStructured& cellSet) +{ + out << "DATASET RECTILINEAR_GRID\n"; + + WriteDimensions(out, cellSet); + + std::string typeName = vtkm::io::internal::DataTypeName::Name(); + vtkm::cont::ArrayHandle dimArray; + + dimArray = points.GetStorage().GetFirstArray(); + out << "X_COORDINATES " << dimArray.GetNumberOfValues() << " " << typeName << "\n"; + OutputFieldFunctor{ out }(dimArray); + + dimArray = points.GetStorage().GetSecondArray(); + out << "Y_COORDINATES " << dimArray.GetNumberOfValues() << " " << typeName << "\n"; + OutputFieldFunctor{ out }(dimArray); + + dimArray = points.GetStorage().GetThirdArray(); + out << "Z_COORDINATES " << dimArray.GetNumberOfValues() << " " << typeName << "\n"; + OutputFieldFunctor{ out }(dimArray); +} + +template +void WriteDataSetAsStructuredGrid(std::ostream& out, + const vtkm::cont::DataSet& dataSet, + const vtkm::cont::CellSetStructured& cellSet) +{ + out << "DATASET STRUCTURED_GRID" << '\n'; + + WriteDimensions(out, cellSet); + + WritePoints(out, dataSet); +} + template void WriteDataSetAsStructured(std::ostream& out, const vtkm::cont::DataSet& dataSet, const vtkm::cont::CellSetStructured& cellSet) { - ///\todo: support uniform/rectilinear - out << "DATASET STRUCTURED_GRID" << '\n'; + ///\todo: support rectilinear - auto pointDimensions = cellSet.GetPointDimensions(); - using VTraits = vtkm::VecTraits; - - out << "DIMENSIONS "; - out << VTraits::GetComponent(pointDimensions, 0) << " "; - out << (DIM > 1 ? VTraits::GetComponent(pointDimensions, 1) : 1) << " "; - out << (DIM > 2 ? VTraits::GetComponent(pointDimensions, 2) : 1) << " "; - - WritePoints(out, dataSet); + // Type of structured grid (uniform, rectilinear, curvilinear) is determined by coordinate system + vtkm::cont::ArrayHandleVirtualCoordinates coordSystem = dataSet.GetCoordinateSystem().GetData(); + if (coordSystem.IsType()) + { + // uniform is written as "structured points" + WriteDataSetAsStructuredPoints( + out, coordSystem.Cast(), cellSet); + } + else if (coordSystem.IsType>()) + { + WriteDataSetAsRectilinearGrid( + out, coordSystem.Cast>(), cellSet); + } + else if (coordSystem.IsType>()) + { + WriteDataSetAsRectilinearGrid( + out, coordSystem.Cast>(), cellSet); + } + else + { + // Curvilinear is written as "structured grid" + WriteDataSetAsStructuredGrid(out, dataSet, cellSet); + } } void Write(std::ostream& out, const vtkm::cont::DataSet& dataSet, bool just_points = false) diff --git a/vtkm/io/VTKStructuredPointsReader.cxx b/vtkm/io/VTKStructuredPointsReader.cxx index 1293da374..442fa37de 100644 --- a/vtkm/io/VTKStructuredPointsReader.cxx +++ b/vtkm/io/VTKStructuredPointsReader.cxx @@ -32,41 +32,66 @@ void VTKStructuredPointsReader::Read() throw vtkm::io::ErrorIO("Incorrect DataSet type"); } - std::string tag; - // Read structured points specific meta-data vtkm::Id3 dim; - vtkm::Vec3f_32 origin, spacing; + vtkm::Vec3f origin; + vtkm::Vec3f spacing; - //Two ways the file can describe the dimensions. The proper way is by - //using the DIMENSIONS keyword, but VisIt written VTK files spicify data - //bounds instead, as a FIELD + // The specification for VTK Legacy files says the order of fields should be + // DIMENSIONS, ORIGIN, SPACING. However, it is common for these to be in + // different orders. In particular, SPACING is often put before ORIGIN (even + // in VTK's writer). Also, VisIt writes the DIMENSIONS in a different way. + // Account for these differences. + + bool readDim = false; + bool readOrigin = false; + bool readSpacing = false; std::vector visitBounds; - this->DataFile->Stream >> tag; - if (tag == "FIELD") - { - this->ReadGlobalFields(&visitBounds); - this->DataFile->Stream >> tag; - } - if (visitBounds.empty()) - { - internal::parseAssert(tag == "DIMENSIONS"); - this->DataFile->Stream >> dim[0] >> dim[1] >> dim[2] >> std::ws; - this->DataFile->Stream >> tag; - } - internal::parseAssert(tag == "SPACING"); - this->DataFile->Stream >> spacing[0] >> spacing[1] >> spacing[2] >> std::ws; - if (!visitBounds.empty()) + while (!(readDim && readOrigin && readSpacing)) { - //now with spacing and physical bounds we can back compute the dimensions - dim[0] = static_cast((visitBounds[1] - visitBounds[0]) / spacing[0]); - dim[1] = static_cast((visitBounds[3] - visitBounds[2]) / spacing[1]); - dim[2] = static_cast((visitBounds[5] - visitBounds[4]) / spacing[2]); - } + std::string tag; + this->DataFile->Stream >> tag; + if (tag == "DIMENSIONS") + { + this->DataFile->Stream >> dim[0] >> dim[1] >> dim[2] >> std::ws; + readDim = true; + } + else if (tag == "ORIGIN") + { + this->DataFile->Stream >> origin[0] >> origin[1] >> origin[2] >> std::ws; + readOrigin = true; + } + else if (tag == "SPACING") + { + this->DataFile->Stream >> spacing[0] >> spacing[1] >> spacing[2] >> std::ws; + readSpacing = true; + } + else if (tag == "FIELD") + { + // VisIt adds its own metadata in a FIELD tag. + this->ReadGlobalFields(&visitBounds); + } + else + { + std::stringstream message("Parse error: unexpected tag "); + message << tag; + throw vtkm::io::ErrorIO(message.str()); + } - this->DataFile->Stream >> tag >> origin[0] >> origin[1] >> origin[2] >> std::ws; - internal::parseAssert(tag == "ORIGIN"); + //Two ways the file can describe the dimensions. The proper way is by + //using the DIMENSIONS keyword, but VisIt written VTK files spicify data + //bounds instead, as a FIELD + if (readSpacing && !visitBounds.empty()) + { + //now with spacing and physical bounds we can back compute the dimensions + dim[0] = static_cast((visitBounds[1] - visitBounds[0]) / spacing[0]); + dim[1] = static_cast((visitBounds[3] - visitBounds[2]) / spacing[1]); + dim[2] = static_cast((visitBounds[5] - visitBounds[4]) / spacing[2]); + readDim = true; + visitBounds.clear(); + } + } this->DataSet.SetCellSet(internal::CreateCellSetStructured(dim)); this->DataSet.AddCoordinateSystem( diff --git a/vtkm/io/testing/UnitTestVTKDataSetWriter.cxx b/vtkm/io/testing/UnitTestVTKDataSetWriter.cxx index 6925b6e3c..1e94e778b 100644 --- a/vtkm/io/testing/UnitTestVTKDataSetWriter.cxx +++ b/vtkm/io/testing/UnitTestVTKDataSetWriter.cxx @@ -11,6 +11,7 @@ #include #include #include +#include #include #include @@ -22,11 +23,107 @@ namespace #define WRITE_FILE(MakeTestDataMethod) \ TestVTKWriteTestData(#MakeTestDataMethod, tds.MakeTestDataMethod()) +struct CheckSameField +{ + template + void operator()(const vtkm::cont::ArrayHandle& originalArray, + const vtkm::cont::Field& fileField) const + { + vtkm::cont::ArrayHandle fileArray; + fileField.GetData().CopyTo(fileArray); + VTKM_TEST_ASSERT(test_equal_portals(originalArray.ReadPortal(), fileArray.ReadPortal())); + } +}; + +struct CheckSameCoordinateSystem +{ + template + void operator()(const vtkm::cont::ArrayHandle& originalArray, + const vtkm::cont::CoordinateSystem& fileCoords) const + { + CheckSameField{}(originalArray, fileCoords); + } + + void operator()(const vtkm::cont::ArrayHandleUniformPointCoordinates& originalArray, + const vtkm::cont::CoordinateSystem& fileCoords) const + { + VTKM_TEST_ASSERT(fileCoords.GetData().IsType()); + vtkm::cont::ArrayHandleUniformPointCoordinates fileArray = + fileCoords.GetData().Cast(); + auto originalPortal = originalArray.ReadPortal(); + auto filePortal = fileArray.ReadPortal(); + VTKM_TEST_ASSERT(test_equal(originalPortal.GetOrigin(), filePortal.GetOrigin())); + VTKM_TEST_ASSERT(test_equal(originalPortal.GetSpacing(), filePortal.GetSpacing())); + VTKM_TEST_ASSERT(test_equal(originalPortal.GetRange3(), filePortal.GetRange3())); + } + + using ArrayHandleRectilinearCoords = + vtkm::cont::ArrayHandleCartesianProduct, + vtkm::cont::ArrayHandle, + vtkm::cont::ArrayHandle>; + void operator()(const ArrayHandleRectilinearCoords& originalArray, + const vtkm::cont::CoordinateSystem& fileCoords) const + { + VTKM_TEST_ASSERT(fileCoords.GetData().IsType()); + ArrayHandleRectilinearCoords fileArray = + fileCoords.GetData().Cast(); + auto originalPortal = originalArray.ReadPortal(); + auto filePortal = fileArray.ReadPortal(); + VTKM_TEST_ASSERT( + test_equal_portals(originalPortal.GetFirstPortal(), filePortal.GetFirstPortal())); + VTKM_TEST_ASSERT( + test_equal_portals(originalPortal.GetSecondPortal(), filePortal.GetSecondPortal())); + VTKM_TEST_ASSERT( + test_equal_portals(originalPortal.GetThirdPortal(), filePortal.GetThirdPortal())); + } + + void operator()(const vtkm::cont::ArrayHandleVirtualCoordinates& originalArray, + const vtkm::cont::CoordinateSystem& fileCoords) const + { + if (originalArray.IsType()) + { + (*this)(originalArray.Cast(), fileCoords); + } + else if (originalArray.IsType()) + { + (*this)(originalArray.Cast(), fileCoords); + } + else + { + CheckSameField{}(originalArray, fileCoords); + } + } +}; + +void CheckWrittenReadData(const vtkm::cont::DataSet& originalData, + const vtkm::cont::DataSet& fileData) +{ + VTKM_TEST_ASSERT(originalData.GetNumberOfPoints() == fileData.GetNumberOfPoints()); + VTKM_TEST_ASSERT(originalData.GetNumberOfCells() == fileData.GetNumberOfCells()); + + for (vtkm::IdComponent fieldId = 0; fieldId < originalData.GetNumberOfFields(); ++fieldId) + { + vtkm::cont::Field originalField = originalData.GetField(fieldId); + VTKM_TEST_ASSERT(fileData.HasField(originalField.GetName(), originalField.GetAssociation())); + vtkm::cont::Field fileField = + fileData.GetField(originalField.GetName(), originalField.GetAssociation()); + vtkm::cont::CastAndCall(originalField, CheckSameField{}, fileField); + } + + VTKM_TEST_ASSERT(fileData.GetNumberOfCoordinateSystems() > 0); + CheckSameCoordinateSystem{}(originalData.GetCoordinateSystem().GetData(), + fileData.GetCoordinateSystem()); +} + void TestVTKWriteTestData(const std::string& methodName, const vtkm::cont::DataSet& data) { std::cout << "Writing " << methodName << std::endl; vtkm::io::VTKDataSetWriter writer(methodName + ".vtk"); writer.WriteDataSet(data); + + // Read back and check. + vtkm::io::VTKDataSetReader reader(methodName + ".vtk"); + CheckWrittenReadData(data, reader.ReadDataSet()); } void TestVTKExplicitWrite()