Writing VTK XML files in C++
I visualize the output of my simulations using either LLNL’s Visit or Kitware’s ParaView. These tools are wonderful for dealing with large datasets and can read a huge variety of file formats. In particular, they are good for remote and visualization.
One of the formats that both these tools can understand is the VTK XML format. The advantage of this format is its simplicity and that makes it suitable for small research codes.
Note: If your code is running on large modern HPC machines, you should consider formats that are more efficient for I/O and remote visualization.
In this blog post I’ll describe how I write data from my particle-based simulations into VTK format output files. If the files are small, I use ASCII for readability but for large simulations I write out VTK XML files containing binary data.
The data
The data produced by the particle simulation codes typically has two components:
- A grid that describes the boundary of the computational domain and the partitioning of the domain into pieces that are sent to various processors during the computation process, and
- A set of particles that contain position, geometry, and physical state data.
Each of these datasets changes with time; so they also contain a timestamp. A set of files is produced at each timestep in the simulation where the user feels that there is a need to output results. The process of creating appropriately named files is routine and won’t be described in this post.
Installing VTK and setting up Cmake
In Ubuntu 16.04, installing VTK (a slightly older version, 5.10), is straightforward. I haven’t tried more recent versions, but the API hasn’t changed much. Most HPC systems also have VTK installed by default because of its wide use. Installation just needs
sudo apt-get install libvtk5-dev
To be able to build your code with VTK support, edit the root CMakeLists.txt
file and add
find_package(VTK 5.10 REQUIRED)
include_directories(${VTK_INCLUDE_DIRS})
and for linking the executable, add
target_link_libraries(your_executable
YOUR_LIB
${VTK_LIBRARIES}
)
Writing grid data
We write the grid using an unstructured grid writer (because the grid may not be regular) but stick to hexahedral elements in this example.
#include <vtkXMLUnstructuredGrid.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <vtkSmartPointer.h>
#include <vtkPoints.h>
using vtkPointsP = vtkSmartPointer<vtkPoints>;
using vtkUnstructuredGridP = vtkSmartPointer<vtkUnstructuredGrid>;
using vtkHexahedronP = vtkSmartPointer<vtkHexahedron>;
using vtkXMLUnstructuredGridWriterP = vtkSmartPointer<vtkXMLUnstructuredGridWriter>;
using Node = std::vector<double>;
using Element = std::vector<int>;
void
OutputVTK::writeGrid(double time,
const BoxArray& grid,
std::ostringstream& fileName)
{
// Create a writer
auto writer = vtkXMLUnstructuredGridWriterP::New();
// Append the default extension to the file name
fileName << "." << writer->GetDefaultFileExtension();
writer->SetFileName((fileName.str()).c_str());
// Create a pointer to a VTK Unstructured Grid data set
auto dataSet = vtkUnstructuredGridP::New();
// Set up pointer to point data
auto pts = vtkPointsP::New();
// Count the total number of points to be saved
int num_pts = grid.getNumberOfPoints(); // Implementation is user-dependent
pts->SetNumberOfPoints(num_pts);
// Add the time
addTimeToVTKDataSet(time, dataSet);
// Get the nodes and elements that are used to describe the grid
std::vector<Node> nodes = grid.getNodes(); // Implementation is user-dependent
std::vector<Element> elements = grid.getElements(); // Implementation is user-dependent
// Add the processor boundaries to the unstructured grid cell data
addElementsToVTKDataSet(nodes, elements, pts, dataSet);
// Set the points
dataSet->SetPoints(pts);
// Remove unused memory
dataSet->Squeeze();
// Write the data
writer->SetInput(dataSet);
writer->SetDataModeToAscii();
writer->Write();
}
We will discuss the implementation of addTimeToVTKDataSet
and addElementsToVTKDataSet
next. Note the SetDataModeToAscii()
function near the end. We can switch to binary output by just changing this call to SetDataModeToBinary()
.
Adding time to the data set
To add a time stamp to each data set, we use the following
void
OutputVTK::addTimeToVTKDataSet(double time,
vtkUnstructuredGridP& dataSet)
{
auto array = vtkDoubleArrayP::New();
array->SetName("TIME");
array->SetNumberOfTuples(1);
array->SetTuple1(0, time);
dataSet->GetFieldData()->AddArray(array);
}
Adding elements to the data set
Finally, we will add the grid elements to the dataset using the addElementsToVTKDataSet
function below.
void
OutputVTK::addElementsToVTKDataSet(const std::vector<Node>& nodes,
const std::vector<Element>& elements,
vtkPointsP& pts,
vtkUnstructuredGridP& dataSet)
{
// Set the coordinates of the nodes
int id = 0;
for (const auto& node : nodes) {
pts->SetPoint(id, node[0], node[1], node[2]);
++id;
}
// Set the element connectivities
auto hex = vtkHexahedronP::New(); // Assuming hex elements
for (const auto& element : elements) {
// Get node ids and assign them to a hex element (ids start from 1)
int nodeNum = 0;
for (const auto& id : element) {
hex->GetPointIds()->SetId(nodeNum, id - 1);
++nodeNum;
}
dataSet->InsertNextCell(hex->GetCellType(), hex->GetPointIds());
}
}
The output file
The output file produced by this approach has the extension .vtu
and the format below.
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian" compressor="vtkZLibDataCompressor">
<UnstructuredGrid>
<FieldData>
<DataArray type="Float64" Name="TIME" NumberOfTuples="1" format="ascii" RangeMin="0" RangeMax="0">
0
</DataArray>
</FieldData>
<Piece NumberOfPoints="18" NumberOfCells="4">
<PointData>
</PointData>
<CellData>
</CellData>
<Points>
<DataArray type="Float32" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="101.56425869">
0 -1 0 0 -1 17.728469849
0 0 0 0 0 17.728469849
0 1 0 0 1 17.728469849
50 -1 0 50 -1 17.728469849
50 0 0 50 0 17.728469849
50 1 0 50 1 17.728469849
100 -1 0 100 -1 17.728469849
100 0 0 100 0 17.728469849
100 1 0 100 1 17.728469849
</DataArray>
</Points>
<Cells>
<DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="17">
6 8 2 0 7 9
3 1 8 10 4 2
9 11 5 3 12 14
8 6 13 15 9 7
14 16 10 8 15 17
11 9
</DataArray>
<DataArray type="Int64" Name="offsets" format="ascii" RangeMin="8" RangeMax="32">
8 16 24 32
</DataArray>
<DataArray type="UInt8" Name="types" format="ascii" RangeMin="12" RangeMax="12">
12 12 12 12
</DataArray>
</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>
The data can be easily read into Visit and the resulting plot can be seen below.
Remarks
Writing your data to files in VTK format is straightforward and reading them in Visit or ParaView is also quite simple. I prefer these tools to TecPlot, mainly because I can write plugins quite easily.
In the next post, I will discuss how you can save your particle data in VTK binary format.
If you have questions/comments/corrections, please contact banerjee at parresianz dot com dot zen (without the dot zen).