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 should be (on 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 patch will fix this problem, but is this an issue of VTK?
diff --git a/IO/Xdmf3/vtkXdmf3DataSet.cxx b/IO/Xdmf3/vtkXdmf3DataSet.cxx
index f595132ade..2d9e323c4e 100644
--- a/IO/Xdmf3/vtkXdmf3DataSet.cxx
+++ b/IO/Xdmf3/vtkXdmf3DataSet.cxx
@@ -1341,7 +1341,7 @@ void vtkXdmf3DataSet::VTKToXdmf(
bool OK = true;
vtkDataArray *vCoords = dataSet->GetXCoordinates();
- OK &= vtkXdmf3DataSet::VTKToXdmfArray(vCoords, xZCoords.get());
+ OK &= vtkXdmf3DataSet::VTKToXdmfArray(vCoords, xXCoords.get());
if (OK)
{
vCoords = dataSet->GetYCoordinates();
@@ -1349,7 +1349,7 @@ void vtkXdmf3DataSet::VTKToXdmf(
if (OK)
{
vCoords = dataSet->GetZCoordinates();
- OK &= vtkXdmf3DataSet::VTKToXdmfArray(vCoords, xXCoords.get());
+ OK &= vtkXdmf3DataSet::VTKToXdmfArray(vCoords, xZCoords.get());
}
}