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).

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).

vtkm/io/VTKDataSetWriter.cxx

@@ -34,6 +34,12 @@
 namespace
 {
 
+template <typename T>
+using ArrayHandleRectilinearCoordinates =
+  vtkm::cont::ArrayHandleCartesianProduct<vtkm::cont::ArrayHandle<T>,
+                                          vtkm::cont::ArrayHandle<T>,
+                                          vtkm::cont::ArrayHandle<T>>;
+
 struct OutputPointsFunctor
 {
 private:
@@ -135,6 +141,18 @@ private:
   std::string* Name;
 };
 
+template <vtkm::IdComponent DIM>
+void WriteDimensions(std::ostream& out, const vtkm::cont::CellSetStructured<DIM>& cellSet)
+{
+  auto pointDimensions = cellSet.GetPointDimensions();
+  using VTraits = vtkm::VecTraits<decltype(pointDimensions)>;
+
+  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 <vtkm::IdComponent DIM>
+void WriteDataSetAsStructuredPoints(std::ostream& out,
+                                    const vtkm::cont::ArrayHandleUniformPointCoordinates& points,
+                                    const vtkm::cont::CellSetStructured<DIM>& 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 <typename T, vtkm::IdComponent DIM>
+void WriteDataSetAsRectilinearGrid(std::ostream& out,
+                                   const ArrayHandleRectilinearCoordinates<T>& points,
+                                   const vtkm::cont::CellSetStructured<DIM>& cellSet)
+{
+  out << "DATASET RECTILINEAR_GRID\n";
+
+  WriteDimensions(out, cellSet);
+
+  std::string typeName = vtkm::io::internal::DataTypeName<T>::Name();
+  vtkm::cont::ArrayHandle<T> 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 <vtkm::IdComponent DIM>
+void WriteDataSetAsStructuredGrid(std::ostream& out,
+                                  const vtkm::cont::DataSet& dataSet,
+                                  const vtkm::cont::CellSetStructured<DIM>& cellSet)
+{
+  out << "DATASET STRUCTURED_GRID" << '\n';
+
+  WriteDimensions(out, cellSet);
+
+  WritePoints(out, dataSet);
+}
+
 template <vtkm::IdComponent DIM>
 void WriteDataSetAsStructured(std::ostream& out,
                               const vtkm::cont::DataSet& dataSet,
                               const vtkm::cont::CellSetStructured<DIM>& cellSet)
 {
-  ///\todo: support uniform/rectilinear
-  out << "DATASET STRUCTURED_GRID" << '\n';
+  ///\todo: support rectilinear
 
-  auto pointDimensions = cellSet.GetPointDimensions();
-  using VTraits = vtkm::VecTraits<decltype(pointDimensions)>;
-
-  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<vtkm::cont::ArrayHandleUniformPointCoordinates>())
+  {
+    // uniform is written as "structured points"
+    WriteDataSetAsStructuredPoints(
+      out, coordSystem.Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>(), cellSet);
+  }
+  else if (coordSystem.IsType<ArrayHandleRectilinearCoordinates<vtkm::Float32>>())
+  {
+    WriteDataSetAsRectilinearGrid(
+      out, coordSystem.Cast<ArrayHandleRectilinearCoordinates<vtkm::Float32>>(), cellSet);
+  }
+  else if (coordSystem.IsType<ArrayHandleRectilinearCoordinates<vtkm::Float64>>())
+  {
+    WriteDataSetAsRectilinearGrid(
+      out, coordSystem.Cast<ArrayHandleRectilinearCoordinates<vtkm::Float64>>(), 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)

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<vtkm::Float32> 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<vtkm::Id>((visitBounds[1] - visitBounds[0]) / spacing[0]);
-    dim[1] = static_cast<vtkm::Id>((visitBounds[3] - visitBounds[2]) / spacing[1]);
-    dim[2] = static_cast<vtkm::Id>((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<vtkm::Id>((visitBounds[1] - visitBounds[0]) / spacing[0]);
+      dim[1] = static_cast<vtkm::Id>((visitBounds[3] - visitBounds[2]) / spacing[1]);
+      dim[2] = static_cast<vtkm::Id>((visitBounds[5] - visitBounds[4]) / spacing[2]);
+      readDim = true;
+      visitBounds.clear();
+    }
+  }
 
   this->DataSet.SetCellSet(internal::CreateCellSetStructured(dim));
   this->DataSet.AddCoordinateSystem(

vtkm/io/testing/UnitTestVTKDataSetWriter.cxx

@@ -11,6 +11,7 @@
 #include <complex>
 #include <cstdio>
 #include <vector>
+#include <vtkm/io/VTKDataSetReader.h>
 #include <vtkm/io/VTKDataSetWriter.h>
 
 #include <vtkm/cont/testing/MakeTestDataSet.h>
@@ -22,11 +23,107 @@ namespace
 #define WRITE_FILE(MakeTestDataMethod)                                                             \
   TestVTKWriteTestData(#MakeTestDataMethod, tds.MakeTestDataMethod())
 
+struct CheckSameField
+{
+  template <typename T, typename S>
+  void operator()(const vtkm::cont::ArrayHandle<T, S>& originalArray,
+                  const vtkm::cont::Field& fileField) const
+  {
+    vtkm::cont::ArrayHandle<T> fileArray;
+    fileField.GetData().CopyTo(fileArray);
+    VTKM_TEST_ASSERT(test_equal_portals(originalArray.ReadPortal(), fileArray.ReadPortal()));
+  }
+};
+
+struct CheckSameCoordinateSystem
+{
+  template <typename T>
+  void operator()(const vtkm::cont::ArrayHandle<T>& 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>());
+    vtkm::cont::ArrayHandleUniformPointCoordinates fileArray =
+      fileCoords.GetData().Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>();
+    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::FloatDefault>,
+                                            vtkm::cont::ArrayHandle<vtkm::FloatDefault>,
+                                            vtkm::cont::ArrayHandle<vtkm::FloatDefault>>;
+  void operator()(const ArrayHandleRectilinearCoords& originalArray,
+                  const vtkm::cont::CoordinateSystem& fileCoords) const
+  {
+    VTKM_TEST_ASSERT(fileCoords.GetData().IsType<ArrayHandleRectilinearCoords>());
+    ArrayHandleRectilinearCoords fileArray =
+      fileCoords.GetData().Cast<ArrayHandleRectilinearCoords>();
+    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<vtkm::cont::ArrayHandleUniformPointCoordinates>())
+    {
+      (*this)(originalArray.Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>(), fileCoords);
+    }
+    else if (originalArray.IsType<ArrayHandleRectilinearCoords>())
+    {
+      (*this)(originalArray.Cast<ArrayHandleRectilinearCoords>(), 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()