Skip to content

Fix assignment of coordinates for XdmfRectilnearGrid

For the following program, which writes a RectilinearGrid dataset with Xdmf3Writer, I found that the coordinate arrays in the geometry type of VXVYVZ are written in order of Z, Y and then X.

#include "vtkSmartPointer.h"
#include "vtkFloatArray.h"
#include "vtkCellData.h"
#include "vtkRectilinearGrid.h"
#include "vtkXdmf3Writer.h"

const static double x[5] = {
  0.0, 0.5, 1.0, 1.5, 2.0
};

const static double y[4] = {
  0.0, 0.3, 0.7, 1.0,
};

const static double z[5] = {
  0.0, 0.2, 0.4, 0.6, 0.8,
};

const static double d[2*3*4] = {
  0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0,
  7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0,
  14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0,
  21.0, 22.0, 23.0,
};

int main(int argc, char **argv)
{
  vtkFloatArray *xCoords = vtkFloatArray::New();
  vtkFloatArray *yCoords = vtkFloatArray::New();
  vtkFloatArray *zCoords = vtkFloatArray::New();
  vtkFloatArray *dArray = vtkFloatArray::New();
  vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
  vtkSmartPointer<vtkXdmf3Writer> writer =
    vtkSmartPointer<vtkXdmf3Writer>::New();

  for (int i = 0; i < 5; ++i) xCoords->InsertNextValue(x[i]);
  for (int i = 0; i < 4; ++i) yCoords->InsertNextValue(y[i]);
  for (int i = 0; i < 3; ++i) zCoords->InsertNextValue(z[i]);
  for (int i = 0; i < 24; ++i) dArray->InsertNextValue(d[i]);

  xCoords->SetName("x coords");
  yCoords->SetName("y coords");
  zCoords->SetName("z coords");
  dArray->SetName("t");

  rGrid->SetDimensions(3, 4, 5);
  rGrid->SetXCoordinates(xCoords);
  rGrid->SetYCoordinates(yCoords);
  rGrid->SetZCoordinates(zCoords);
  rGrid->GetCellData()->AddArray(dArray);

  writer->SetFileName("test-single.xmf");
  writer->SetInputData(rGrid);
  writer->Update();
  writer->Write();

  rGrid->Delete();
  xCoords->Delete();
  yCoords->Delete();
  zCoords->Delete();
  dArray->Delete();

  return 0;
}

The output will be (in VTK 8.1.1):

<?xml version="1.0" encoding="utf-8"?>
<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="3.0">
  <Domain>
    <Grid Name="Grid">
      <Geometry Origin="" Type="VXVYVZ">
        <DataItem DataType="Float" Dimensions="5" Format="XML" Name="z coords" Precision="4">0 0.2 0.40000001 0.60000002 0.80000001</DataItem>
        <DataItem DataType="Float" Dimensions="4" Format="XML" Name="y coords" Precision="4">0 0.30000001 0.69999999 1</DataItem>
        <DataItem DataType="Float" Dimensions="3" Format="XML" Name="x coords" Precision="4">0 0.5 1</DataItem>
      </Geometry>
      <Topology Dimensions="5 4 3" Type="3DRectMesh"/>
      <Attribute Center="Cell" ElementCell="" ElementDegree="0" ElementFamily="" ItemType="" Name="t" Type="None">
        <DataItem DataType="Float" Dimensions="4 3 2" Format="XML" Precision="4">0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23</DataItem>
      </Attribute>
    </Grid>
  </Domain>
</Xdmf>

ParaView (at least) assigns the first item to X axis, the second item to Y axis, and the third item to Z axis. So this will be a different result to the dateset written by XMLRectilinearGridWriter.

I found that VTK passes the XCoodinates as ZCoords and ZCoordinates as XCoords for XDMF library, so this change set will fix this problem.

This report was originally posted to discourse.

PS. I wrote this commit message, but I could not find any source about this (The XDMF specification says that the indices of I varies at first, but it does not seem to say that it applies to X coordinates).

Both XdmfRectilinearGrid and vtkRectilinearGrid use the topology that
X (or I) indices varies at first.

Merge request reports