Cookbook

Short answers to "how do I ...?". Each recipe is a complete program with its real output; the reference has the details, the tutorial the whole story.

• Write each kind of dataset
  • A regular grid (ImageData)
  • A curvilinear grid (StructuredGrid)
  • Cells of any type, polyhedra included (UnstructuredGrid)
  • Points, lines and polygons (PolyData)
• Write data
  • Scalars, vectors and tensors
  • Time, cycle and other global data
  • Ghost cells and other unsigned arrays
• Choose the format
  • Compress, and go beyond 2 GiB
  • Write a .vtkhdf file
  • A PolyData in a .vtkhdf file
  • A grid in a .vtkhdf file
  • A time series in one .vtkhdf file
  • A fixed mesh stored once
  • A time series of several meshes in one .vtkhdf file
  • A time history (Table)
  • A block-structured AMR
  • A tree-based AMR
  • More than 2^31 points or cells
  • Several pieces in one file
• Parallel and composite files
  • A partitioned structured grid (.pvts)
  • Write from MPI ranks
  • Write one .vtkhdf from MPI ranks
  • Check a parallel header against its pieces
  • Continue a time series after a restart
  • Assemblies of datasets (.vtm)
  • Write a file into memory
• Read files
  • What a file holds
  • Read one array
  • Read the mesh
  • Handle the errors

Write each kind of dataset

A regular grid (ImageData)

No coordinates at all: the points are defined by the extent, an origin and a spacing.

program image_data
!< Write a regular grid (ImageData): no coordinates, just an origin and a spacing.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
integer(I4P), parameter :: nx=40, ny=30, nz=10
real(R8P)               :: wave(0:nx,0:ny,0:nz)
type(vtk_file)          :: a_vtk_file
integer(I4P)            :: i, j, k, error

do k=0, nz ; do j=0, ny ; do i=0, nx
  wave(i,j,k) = sin(0.3_R8P*i)*cos(0.4_R8P*j) + 0.1_R8P*k
enddo ; enddo ; enddo
error = a_vtk_file%initialize(format='raw', filename='wave.vti', mesh_topology='ImageData', &
                              nx1=0, nx2=nx, ny1=0, ny2=ny, nz1=0, nz2=nz,                  &
                              origin=[0._R8P, 0._R8P, 0._R8P], spacing=[0.1_R8P, 0.1_R8P, 0.1_R8P])
error = a_vtk_file%xml_writer%write_piece(nx1=0, nx2=nx, ny1=0, ny2=ny, nz1=0, nz2=nz)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='wave', x=wave, one_component=.true.)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
print '(A,I0)', 'wave.vti written, error ', error
endprogram image_data

$ image_data
wave.vti written, error 0

direction (a row-major 3x3 matrix) rotates the axes of the grid. In a parallel .pvti, every piece has the same origin and spacing.

A curvilinear grid (StructuredGrid)

Every point has its own coordinates, the topology is still that of a box: write_geo(n, x, y, z) with rank-3 arrays.

program curvilinear
!< Write a curvilinear grid (StructuredGrid): every point has its own coordinates, here a quarter of an annulus.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
integer(I4P), parameter :: nr=8, nt=24, nz=2
real(R8P),    parameter :: pi=acos(-1._R8P)
real(R8P)               :: x(nr,nt,nz), y(nr,nt,nz), z(nr,nt,nz), r(nr,nt,nz)
type(vtk_file)          :: a_vtk_file
integer(I4P)            :: i, j, k, error

do k=1, nz ; do j=1, nt ; do i=1, nr
  r(i,j,k) = 1 + (i - 1)/real(nr - 1, R8P)
  x(i,j,k) = r(i,j,k)*cos(pi/2*(j - 1)/(nt - 1))
  y(i,j,k) = r(i,j,k)*sin(pi/2*(j - 1)/(nt - 1))
  z(i,j,k) = 0.1_R8P*(k - 1)
enddo ; enddo ; enddo
error = a_vtk_file%initialize(format='binary', filename='annulus.vts', mesh_topology='StructuredGrid', &
                              nx1=1, nx2=nr, ny1=1, ny2=nt, nz1=1, nz2=nz)
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=nr, ny1=1, ny2=nt, nz1=1, nz2=nz)
error = a_vtk_file%xml_writer%write_geo(n=nr*nt*nz, x=x, y=y, z=z)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='radius', x=r, one_component=.true.)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
print '(A,I0)', 'annulus.vts written, error ', error
endprogram curvilinear

$ curvilinear
annulus.vts written, error 0

Cells of any type, polyhedra included (UnstructuredGrid)

cell_type gives the VTK type of each cell, offset the end of each cell in connect. A polyhedron (type 42) also needs its faces: face and faceoffset, -1 for the cells that are not polyhedra; the points of each face are listed counter-clockwise seen from outside the cell.

program polyhedron
!< Write an unstructured grid of different cells: a polyhedron (a cube described by its faces), a tetrahedron, a wedge.
use penf, only : I1P, I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
real(R8P),    parameter :: x(18)=[real(R8P) :: 0,1,1,0,0,1,1,0, 1.5,2.5,1.5,1.5, 3,4,3,3,4,3]
real(R8P),    parameter :: y(18)=[real(R8P) :: 0,0,1,1,0,0,1,1, 0,0,1,0, 0,0,1,0,0,1]
real(R8P),    parameter :: z(18)=[real(R8P) :: 0,0,0,0,1,1,1,1, 0,0,0,1, 0,0,0,1,1,1]
integer(I4P), parameter :: connect(18)=[0,1,2,3,4,5,6,7, 8,9,10,11, 12,13,14,15,16,17]
integer(I4P), parameter :: offset(3)=[8, 12, 18]
integer(I1P), parameter :: cell_type(3)=[42_I1P, 10_I1P, 13_I1P] ! polyhedron, tetrahedron, wedge
! the faces of the polyhedron: their number, then for each face its number of points and their ids, counter-clockwise
! seen from outside the cell
integer(I4P), parameter :: face(31)=[6, 4,0,3,2,1, 4,4,5,6,7, 4,0,1,5,4, 4,1,2,6,5, 4,2,3,7,6, 4,3,0,4,7]
integer(I4P), parameter :: faceoffset(3)=[31, -1, -1] ! the end of the faces of each cell, -1 if not a polyhedron
type(vtk_file)          :: a_vtk_file
integer(I4P)            :: error

error = a_vtk_file%initialize(format='ascii', filename='cells.vtu', mesh_topology='UnstructuredGrid')
error = a_vtk_file%xml_writer%write_piece(np=18, nc=3)
error = a_vtk_file%xml_writer%write_geo(np=18, nc=3, x=x, y=y, z=z)
error = a_vtk_file%xml_writer%write_connectivity(nc=3, connectivity=connect, offset=offset, cell_type=cell_type, &
                                                 face=face, faceoffset=faceoffset)
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='cell_type', x=int(cell_type, I4P))
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
print '(A,I0)', 'cells.vtu written, error ', error
endprogram polyhedron

$ polyhedron
cells.vtu written, error 0

Points, lines and polygons (PolyData)

Up to four blocks of cells: vertices, lines (polylines), triangle strips, polygons. Pass only the blocks you have.

program polydata
!< Write polygonal data: vertices, a polyline and two polygons, with one value per cell.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
real(R8P)      :: x(12), y(12), z(12)
type(vtk_file) :: a_vtk_file
integer(I4P)   :: error

x = [real(R8P) :: 0, 1, 2,   0, 1, 2, 3,   0, 1, 1,   2, 3]
y = [real(R8P) :: 3, 3, 3,   2, 2.5, 2, 2.5,   0, 0, 1,   0, 1]
z = 0
error = a_vtk_file%initialize(format='raw', filename='shapes.vtp', mesh_topology='PolyData')
error = a_vtk_file%xml_writer%write_piece(np=12, nverts=3, nlines=1, nstrips=0, npolys=2)
error = a_vtk_file%xml_writer%write_geo(np=12, nc=6, x=x, y=y, z=z)
error = a_vtk_file%xml_writer%write_polydata_cells(verts_connectivity=[0, 1, 2], verts_offset=[1, 2, 3],   &
                                                   lines_connectivity=[3, 4, 5, 6], lines_offset=[4],      &
                                                   polys_connectivity=[7, 8, 9, 8, 10, 11, 9], polys_offset=[3, 7])
! one value per cell, in the order vertices, lines, strips, polygons
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='id', x=[1, 1, 1, 2, 3, 4])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
print '(A,I0)', 'shapes.vtp written, error ', error
endprogram polydata

$ polydata
shapes.vtp written, error 0

The cell data follow the VTK order of the blocks: vertices, lines, strips, polygons.

Write data

Scalars, vectors and tensors

A vector by its components, x, y, z; any number of components as a rank-2 array (components, points): 6 for a symmetric tensor (xx, yy, zz, xy, yz, xz), 9 for a full one.

program vectors_tensors
!< Write scalars, vectors and symmetric tensors at the points of a tetrahedron.
use penf, only : I1P, I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
real(R8P), parameter :: x(4)=[0, 1, 0, 0], y(4)=[0, 0, 1, 0], z(4)=[0, 0, 0, 1]
real(R8P)            :: p(4), u(4), v(4), w(4), s(6,4)
type(vtk_file)       :: a_vtk_file
integer(I4P)         :: error

p = [1, 2, 3, 4]                         ! a scalar
u = 1 ; v = x ; w = -y                   ! a vector, by components
s(:,1) = [1, 2, 3, 0, 0, 0] ; s(:,2) = 2*s(:,1) ; s(:,3) = 3*s(:,1) ; s(:,4) = 4*s(:,1) ! xx yy zz xy yz xz
error = a_vtk_file%initialize(format='raw', filename='fields.vtu', mesh_topology='UnstructuredGrid')
error = a_vtk_file%xml_writer%write_piece(np=4, nc=1)
error = a_vtk_file%xml_writer%write_geo(np=4, nc=1, x=x, y=y, z=z)
error = a_vtk_file%xml_writer%write_connectivity(nc=1, connectivity=[0, 1, 2, 3], offset=[4], cell_type=[10_I1P])
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open', scalars='pressure', vectors='velocity', &
                                              tensors='stress')
error = a_vtk_file%xml_writer%write_dataarray(data_name='pressure', x=p)
error = a_vtk_file%xml_writer%write_dataarray(data_name='velocity', x=u, y=v, z=w)
error = a_vtk_file%xml_writer%write_dataarray(data_name='stress', x=s)   ! rank 2: (components, points)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
endprogram vectors_tensors

$ inspect fields.vtu
fields.vtu: UnstructuredGrid, pieces: 1, compressor: none
  piece 1: 4 points, 1 cells
     node pressure         Float64 1 comp., range  1.0000E+00  4.0000E+00
     node velocity         Float64 3 comp., range -1.0000E+00  1.0000E+00
     node stress           Float64 6 comp., range  0.0000E+00  1.2000E+01

scalars, vectors and tensors choose the arrays readers use by default. inspect is the program of the quick start.

Time, cycle and other global data

program field_data
!< Attach global data to a dataset: numbers, arrays, strings; then read them back.
use penf, only : I4P, I8P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file)                :: a_vtk_file
real(R8P),        allocatable :: residuals(:)
character(len=:), allocatable :: species(:), names(:)
integer(I4P)                  :: error

error = a_vtk_file%initialize(format='binary', filename='case.vtr', mesh_topology='RectilinearGrid', &
                              nx1=1, nx2=2, ny1=1, ny2=2, nz1=1, nz2=2)
error = a_vtk_file%xml_writer%write_fielddata(action='open')
error = a_vtk_file%xml_writer%write_fielddata(data_name='TIME', x=0.25_R8P)
error = a_vtk_file%xml_writer%write_fielddata(data_name='CYCLE', x=100_I8P)
error = a_vtk_file%xml_writer%write_fielddata(data_name='residuals', x=[1.e-2_R8P, 1.e-4_R8P, 1.e-6_R8P])
error = a_vtk_file%xml_writer%write_fielddata(data_name='solver', x='my solver v2.1')
error = a_vtk_file%xml_writer%write_fielddata(data_name='species', x=['N2 ', 'O2 ', 'CO2'])
error = a_vtk_file%xml_writer%write_fielddata(action='close')
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=2, ny1=1, ny2=2, nz1=1, nz2=2)
error = a_vtk_file%xml_writer%write_geo(x=[0._R8P, 1._R8P], y=[0._R8P, 1._R8P], z=[0._R8P, 1._R8P])
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

error = a_vtk_file%initialize(filename='case.vtr', action='read')
error = a_vtk_file%xml_reader%get_dataarray_names(location='field', names=names)
print '(A,*(1X,A))', 'field data:', (trim(names(error)), error=1, size(names))
error = a_vtk_file%xml_reader%read_dataarray(location='field', data_name='residuals', x=residuals)
error = a_vtk_file%xml_reader%read_dataarray(location='field', data_name='species', x=species)
print '(A,*(1X,ES8.1))', 'residuals:', residuals
print '(A,*(1X,A))', 'species:', species
error = a_vtk_file%finalize()
endprogram field_data

$ field_data
field data: TIME CYCLE residuals solver species
residuals:  1.0E-02  1.0E-04  1.0E-06
species: N2  O2  CO2

Field data are written right after initialize, before the first piece: scalars and rank-1 arrays of any kind, strings and arrays of strings (trailing blanks trimmed).

Ghost cells and other unsigned arrays

program ghost_cells
!< Mark ghost cells, the duplicates a piece keeps of its neighbours, so that readers do not draw them twice.
use penf, only : I1P, I2P, I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file)            :: a_vtk_file
integer(I1P), allocatable :: flags(:)
integer(I2P), allocatable :: values(:)
integer(I4P)              :: error

error = a_vtk_file%initialize(format='raw', filename='piece.vtr', mesh_topology='RectilinearGrid', &
                              nx1=1, nx2=5, ny1=1, ny2=2, nz1=1, nz2=2)
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=5, ny1=1, ny2=2, nz1=1, nz2=2)
error = a_vtk_file%xml_writer%write_geo(x=[0._R8P, 1._R8P, 2._R8P, 3._R8P, 4._R8P], y=[0._R8P, 1._R8P], &
                                        z=[0._R8P, 1._R8P])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
! vtkGhostType is an unsigned byte: 0 an owned cell, 1 a duplicate (ghost) cell; here the last cell is a ghost
error = a_vtk_file%xml_writer%write_dataarray_unsigned(data_name='vtkGhostType', x=[0_I1P, 0_I1P, 0_I1P, 1_I1P])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

! unsigned arrays read back into the signed kind of the same width (same bits) or into a wider one (same values)
error = a_vtk_file%initialize(filename='piece.vtr', action='read')
error = a_vtk_file%xml_reader%read_dataarray(location='cell', data_name='vtkGhostType', x=flags)
error = a_vtk_file%xml_reader%read_dataarray(location='cell', data_name='vtkGhostType', x=values)
print '(A,4I4,A,4I4)', 'vtkGhostType as I1P:', flags, ', as I2P:', values
error = a_vtk_file%finalize()
endprogram ghost_cells

$ ghost_cells
vtkGhostType as I1P:   0   0   0   1, as I2P:   0   0   0   1

write_dataarray_unsigned writes I1P...I8P arrays as UInt8...UInt64. In vtkGhostType, 1 marks a duplicate cell and 32 a hidden one: ParaView does not draw them. For partitions (one piece per rank), the duplicated boundary points are ghosts and global ids (global_ids= when opening the point or cell data) tell VTK which points are the same: Partitions, ghosts and global ids.

Choose the format

Compress, and go beyond 2 GiB

program compress
!< Compress the binary data with zlib, and use 64-bit headers for arrays beyond 2 GiB; the reader finds both by itself.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
real(R8P)                     :: v(1000)
type(vtk_file)                :: a_vtk_file
character(len=:), allocatable :: header_type, compressor
integer(I4P)                  :: i, error

v = [(sin(i/50._R8P), i=1, size(v))]
error = a_vtk_file%initialize(format='raw', filename='big.vtr', mesh_topology='RectilinearGrid', &
                              nx1=1, nx2=1000, ny1=1, ny2=1, nz1=1, nz2=1,                       &
                              compressor='zlib', header_type='UInt64')
if (error /= 0) error stop 'zlib is not available: build the library with VTKFORTRAN_USE_ZLIB'
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=1000, ny1=1, ny2=1, nz1=1, nz2=1)
error = a_vtk_file%xml_writer%write_geo(x=[(real(i, R8P), i=1, 1000)], y=[0._R8P], z=[0._R8P])
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='v', x=v)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

error = a_vtk_file%initialize(filename='big.vtr', action='read')
error = a_vtk_file%xml_reader%get_info(header_type=header_type, compressor=compressor)
print '(A,A,A,A)', 'big.vtr: header ', header_type, ', compressor ', compressor
error = a_vtk_file%finalize()
endprogram compress

$ compress
big.vtr: header UInt64, compressor zlib

compressor='zlib' needs the library built with zlib: initialize returns an error otherwise, the time to fall back to an uncompressed format. header_type='UInt64' is needed only for arrays of more than 2 GiB. The sizes of each format are measured in chapter 2 of the tutorial.

Write a .vtkhdf file

format='vtkhdf' writes one HDF5 file in the VTKHDF layout, with the calls of the XML formats: ParaView reads it, and so does any HDF5 tool. Here three blocks of hexahedra, each one a piece of an UnstructuredGrid:

error = a_vtk_file%initialize(format='vtkhdf', filename='blocks.vtkhdf', mesh_topology='UnstructuredGrid')
do b=1, 3
  call block_mesh(b, x, y, z, connect, offset, cell_type) ! the hexahedra of block b, point ids from 0 in the block
  error = a_vtk_file%xml_writer%write_piece(np=np, nc=nc)
  error = a_vtk_file%xml_writer%write_geo(np=np, nc=nc, x=x, y=y, z=z)
  error = a_vtk_file%xml_writer%write_connectivity(nc=nc, connectivity=connect, offset=offset, cell_type=cell_type)
  error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
  error = a_vtk_file%xml_writer%write_dataarray(data_name='block', x=[(b, c=1, nc)])
  error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
  error = a_vtk_file%xml_writer%write_piece()
enddo
error = a_vtk_file%finalize()

Each piece is a partition of the file: its counts are one entry of NumberOfPoints, NumberOfCells and NumberOfConnectivityIds, its points, cells and arrays follow those of the previous pieces in datasets shared by all. Point ids count from 0 in each piece, and the offsets start from 0 in each piece, as VTKHDF wants (VTKFortran's offsets are the ends of the cells: the writer converts them). The program reads the counts back with the HDF5 bindings of the library:

$ vtkhdf_pieces
blocks.vtkhdf, one partition per piece:
  NumberOfPoints          18 18 18
  NumberOfCells           4 4 4
  NumberOfConnectivityIds 32 32 32
  Offsets                 0 8 16 24 32 0 8 16 24 32 0 8 16 24 32

ParaView reads the partitions as one partitioned dataset. Every piece must hold the same arrays: an array missing from a piece, or new after the first one, is refused. It needs the library built with HDF5 (initialize returns an error otherwise). Every topology is written, polyhedra included; the layout of the file: VTKHDF.

A PolyData in a .vtkhdf file

The calls of a PolyData, with format='vtkhdf': here a hexagon with a vertex on each corner, then a ribbon (a triangle strip) with a polyline above it, each one a piece.

error = a_vtk_file%initialize(format='vtkhdf', filename='shapes.vtkhdf', mesh_topology='PolyData')
! piece 1: a hexagon and a vertex on each corner
error = a_vtk_file%xml_writer%write_piece(np=6, nverts=6, nlines=0, nstrips=0, npolys=1)
error = a_vtk_file%xml_writer%write_geo(np=6, nc=7, x=cos(angle), y=sin(angle), z=[(0._R8P, i=1, 6)])
error = a_vtk_file%xml_writer%write_polydata_cells(verts_connectivity=[(i, i=0, 5)], verts_offset=[(i, i=1, 6)], &
                                                   polys_connectivity=[(i, i=0, 5)], polys_offset=[6])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='block', x=[(1, i=1, 7)])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
! piece 2: a ribbon (a triangle strip of 10 points) and a polyline above it (5 points)
error = a_vtk_file%xml_writer%write_piece(np=15, nverts=0, nlines=1, nstrips=1, npolys=0)
error = a_vtk_file%xml_writer%write_geo(np=15, nc=2, x=[(2._R8P + i / 2, i=0, 9), (2._R8P + i, i=0, 4)],          &
                                        y=[(-0.4_R8P + 0.8_R8P * mod(i, 2), i=0, 9), (1.3_R8P, i=0, 4)], &
                                        z=[(0._R8P, i=1, 15)])
error = a_vtk_file%xml_writer%write_polydata_cells(lines_connectivity=[(i, i=10, 14)], lines_offset=[5], &
                                                   strips_connectivity=[(i, i=0, 9)], strips_offset=[10])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='block', x=[2, 2])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

Each piece is a partition, and each block of cells a group of the file, with one entry per piece: a block a piece does not have gets zero cells there. The program reads the counts back with the HDF5 bindings of the library:

$ vtkhdf_polydata
shapes.vtkhdf, the cells of each partition by group:
  NumberOfPoints          6 15
  Vertices/NumberOfCells  6  0
  Lines/NumberOfCells     0  1
  Strips/NumberOfCells    0  1
  Polygons/NumberOfCells  1  0

The cell data follow the order of the blocks (vertices, lines, strips, polygons), as in the XML formats.

A grid in a .vtkhdf file

A StructuredGrid (here a thick quarter of an annulus) takes the calls of the XML formats too, in one piece whose extent is the whole extent:

error = a_vtk_file%initialize(format='vtkhdf', filename='annulus.vtkhdf', mesh_topology='StructuredGrid', &
                              nx1=0, nx2=ni, ny1=0, ny2=nj, nz1=0, nz2=nk)
error = a_vtk_file%xml_writer%write_piece(nx1=0, nx2=ni, ny1=0, ny2=nj, nz1=0, nz2=nk)
error = a_vtk_file%xml_writer%write_geo(xyz=xyz)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='radius', x=radius, one_component=.true.)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

The points are one dataset, Points, whatever the layout given to write_geo; the attribute Dimensions holds the points along each axis. The program prints the layout, read back with the HDF5 bindings of the library:

$ vtkhdf_grid
annulus.vtkhdf, shapes in C order (as h5dump and h5py show them):
/VTKHDF                  Type       StructuredGrid
                         Version    2.7
                         Dimensions 13 7 3
/VTKHDF/Points           shape [3,7,13,3]
/VTKHDF/PointData/radius shape [3,7,13]

A RectilinearGrid is written the same way, with write_geo(x=x, y=y, z=z): its coordinates are the datasets XCoordinates, YCoordinates and ZCoordinates. Both are VTKHDF 2.7: VTK 9.7 reads them, ParaView 6.1 does not (its VTK 9.6 knows neither type); use the XML formats (.vtr, .vts) for a ParaView of VTK 9.6 or older.

A time series in one .vtkhdf file

Instead of one file per output and a .pvd collection, the whole series in one .vtkhdf file: each output starts with write_step(time), then its field data and its pieces, written as usual. Here a wave on a 2D grid, written in two runs: the second one reopens the file with action='append', as a restarted simulation would.

! first run: the file is created, steps 0 to 3
error = a_vtk_file%initialize(format='vtkhdf', filename='wave.vtkhdf', mesh_topology='ImageData', &
                              nx1=0, nx2=nx, ny1=0, ny2=ny, nz1=0, nz2=0,                       &
                              origin=[0._R8P, 0._R8P, 0._R8P], spacing=[1._R8P, 1._R8P, 1._R8P])
do s=0, 3
  error = output(time=0.25_R8P * s)
enddo
error = a_vtk_file%finalize()
! second run (a restart): the file is reopened, steps 4 to 7 follow
error = a_vtk_file%initialize(format='vtkhdf', filename='wave.vtkhdf', action='append')
do s=4, 7
  error = output(time=0.25_R8P * s)
enddo
error = a_vtk_file%finalize()

Each output:

function output(time) result(error)
!< One output: the step, its field data, its piece with the wave at that time.
real(R8P), intent(in) :: time  !< Time.
integer(I4P)          :: error !< Error status.
real(R8P)             :: wave(0:nx,0:ny)
integer(I4P)          :: i, j

do j=0, ny ; do i=0, nx
  wave(i,j) = sin(2._R8P * pi * (i / real(nx, R8P) - 0.5_R8P * time)) * cos(pi * (j / real(ny, R8P) - 0.5_R8P))
enddo ; enddo
error = a_vtk_file%xml_writer%write_step(time=time)
error = a_vtk_file%xml_writer%write_fielddata(action='open')
error = a_vtk_file%xml_writer%write_fielddata(data_name='TIME', x=time)
error = a_vtk_file%xml_writer%write_fielddata(action='close')
error = a_vtk_file%xml_writer%write_piece(nx1=0, nx2=nx, ny1=0, ny2=ny, nz1=0, nz2=0)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='wave', x=wave, one_component=.true.)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
endfunction output

$ vtkhdf_series
wave.vtkhdf: 8 steps, at the times
 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75

• The file is valid after every piece: a run killed in the middle leaves the steps written so far, as a .pvd collection does. action='append' continues the file, and drops what an interrupted run left after its last complete piece. The times must increase: a restart from an earlier time is refused.
• The first step writes every array. A grid (ImageData, RectilinearGrid, StructuredGrid) has one piece per step; the parts of an UnstructuredGrid or a PolyData can change from step to step, and so can the mesh. What does not change is stored once: see the next recipe.
• The layout of the file: VTKHDF; the two kinds of series compared: Time series.

A fixed mesh stored once

A simulation on a fixed mesh need not write it at every output: after the first step, a step without pieces reuses the mesh of the step before, and writes its arrays only, in data groups opened directly after write_step. Here a plate of quadrilaterals (an UnstructuredGrid) and its temperature, six steps:

error = a_vtk_file%initialize(format='vtkhdf', filename='plate.vtkhdf', mesh_topology='UnstructuredGrid')
do s=0, n_steps - 1
  error = a_vtk_file%xml_writer%write_step(time=0.25_R8P * s)
  if (s == 0) then ! the first step: the mesh and the temperature
    error = a_vtk_file%xml_writer%write_piece(np=np, nc=nc)
    error = a_vtk_file%xml_writer%write_geo(np=np, nc=nc, x=x, y=y, z=z)
    error = a_vtk_file%xml_writer%write_connectivity(nc=nc, connectivity=connectivity, offset=offset, cell_type=cell_type)
    error = temperature(time=0._R8P)
    error = a_vtk_file%xml_writer%write_piece()
  else             ! the next steps: no piece, the temperature only (the mesh of the step before)
    error = temperature(time=0.25_R8P * s)
  endif
enddo
error = a_vtk_file%finalize()

function temperature(time) result(error)
!< The temperature of the plate at a time: a hot spot moving along x, as point data.
real(R8P), intent(in) :: time  !< Time.
integer(I4P)          :: error !< Error status.

error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='temperature', &
                                              x=exp(-((x - 2 - 8 * time)**2 + (y - ny / 2._R8P)**2) / 8))
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
endfunction temperature

The program writes the same series with the mesh at every step too, and prints what each file stores:

$ vtkhdf_fixed
plate.vtkhdf:   153 points stored,  3672 bytes of coordinates
  PointOffsets of the steps:   0   0   0   0   0   0
every_step.vtkhdf:   918 points stored, 22032 bytes of coordinates
  PointOffsets of the steps:   0 153 306 459 612 765

• The offsets of Steps repeat those of the step before for what a step does not write: here PointOffsets is 0 at every step. The arrays of a step without pieces hold the points (cells) of all the parts of the mesh, in their order.
• An array not written in a step reuses its values of the step before: in a step of a grid, or in a step without pieces; a step with pieces writes every point and cell array. A field data array not written reuses its values too.
• A grid: a RectilinearGrid or StructuredGrid piece without write_geo reuses the coordinates or points of the step before; a step of a grid can also have no piece at all (every array reused).
• action='append' keeps reusing after a restart. The rules: One VTKHDF file instead of a collection.

A time series of several meshes in one .vtkhdf file

Several meshes moving together, in one file: each one is written as its own time series, with the same times, then vtm_file with format='vtkhdf' copies them as the blocks of one multi-block file. Here a wave (a PolyData surface) and a buoy riding it (an UnstructuredGrid of two hexahedra), six steps.

! each mesh is a time series of its own, with the same times
error = a_vtk_file%initialize(format='vtkhdf', filename='sea.vtkhdf', mesh_topology='PolyData')
do s=0, n_steps - 1
  error = sea(time=0.25_R8P * s)
enddo
error = a_vtk_file%finalize()
error = a_vtk_file%initialize(format='vtkhdf', filename='buoy.vtkhdf', mesh_topology='UnstructuredGrid')
do s=0, n_steps - 1
  error = buoy(time=0.25_R8P * s)
enddo
error = a_vtk_file%finalize()
! the multi-block file: the two series copied as its blocks
error = a_vtm_file%initialize(filename='harbour.vtkhdf', format='vtkhdf')
error = a_vtm_file%write_block(filenames='sea.vtkhdf', name='water')
error = a_vtm_file%write_block(filenames='buoy.vtkhdf', name='bodies')
error = a_vtm_file%finalize()

The steps of each mesh, written as for a single series:

function sea(time) result(error)
!< One step of the sea: a surface of nx x ny quadrilaterals, its height a wave travelling along x.
real(R8P), intent(in) :: time  !< Time.
integer(I4P)          :: error !< Error status.
real(R8P)             :: x(0:nx,0:ny), y(0:nx,0:ny), z(0:nx,0:ny)
integer(I4P)          :: i, j, c(4 * nx * ny)

do j=0, ny ; do i=0, nx
  x(i,j) = i ; y(i,j) = j ; z(i,j) = wave(real(i, R8P), time)
enddo ; enddo
do j=0, ny - 1 ; do i=0, nx - 1 ! the corners of each quadrilateral, counter-clockwise
  c(4 * (j * nx + i) + 1:4 * (j * nx + i) + 4) = [i + j * (nx + 1), i + 1 + j * (nx + 1), i + 1 + (j + 1) * (nx + 1), &
                                                  i + (j + 1) * (nx + 1)]
enddo ; enddo
error = a_vtk_file%xml_writer%write_step(time=time)
error = a_vtk_file%xml_writer%write_piece(np=(nx + 1) * (ny + 1), nverts=0, nlines=0, nstrips=0, npolys=nx * ny)
error = a_vtk_file%xml_writer%write_geo(np=(nx + 1) * (ny + 1), nc=nx * ny, x=reshape(x, [size(x)]), &
                                        y=reshape(y, [size(y)]), z=reshape(z, [size(z)]))
error = a_vtk_file%xml_writer%write_polydata_cells(polys_connectivity=c, polys_offset=[(4 * i, i=1, nx * ny)])
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='height', x=reshape(z, [size(z)]))
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
endfunction sea

function buoy(time) result(error)
!< One step of the buoy: two hexahedra (a 2 x 2 x 0.4 box) at x = 6...8, lifted by the wave at its centre.
real(R8P), intent(in) :: time  !< Time.
integer(I4P)          :: error !< Error status.
real(R8P)             :: x(12), y(12), z(12)
integer(I4P)          :: p

do p=0, 11 ! x = 6 + mod(p, 3), y = 3 + 2 mod(p / 3, 2), z = the wave at x = 7, -0.2 or +0.2
  x(p + 1) = 6 + mod(p, 3) ; y(p + 1) = 3 + 2 * mod(p / 3, 2) ; z(p + 1) = wave(7._R8P, time) - 0.2_R8P + 0.4_R8P * (p / 6)
enddo
error = a_vtk_file%xml_writer%write_step(time=time)
error = a_vtk_file%xml_writer%write_piece(np=12, nc=2)
error = a_vtk_file%xml_writer%write_geo(np=12, nc=2, x=x, y=y, z=z)
error = a_vtk_file%xml_writer%write_connectivity(nc=2, connectivity=[0, 1, 4, 3, 6, 7, 10, 9, 1, 2, 5, 4, 7, 8, 11, 10], &
                                                 offset=[8, 16], cell_type=[12_I1P, 12_I1P])
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='height', x=z)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_piece()
endfunction buoy

The file is read back with the library: vtm_file lists its blocks and datasets, and each dataset (its group in the file, the file of the entry) is read by vtk_file with dataset= as a .vtkhdf of its own:

! the entries of the multi-block file, then each dataset (its path in the file) read as a .vtkhdf of its own
error = a_vtm_file%initialize(filename='harbour.vtkhdf', action='read')
error = a_vtm_file%get_entries(level=level, kind=kind, name=name, file=file)
error = a_vtm_file%finalize()
print '(A)', 'harbour.vtkhdf, its blocks and datasets:'
do e=1, size(level)
  if (trim(kind(e)) == 'block') then
    print '(A)', repeat('  ', level(e))//'block '//trim(name(e))
  else
    error = a_vtk_file%initialize(filename='harbour.vtkhdf', action='read', dataset=trim(file(e)))
    error = a_vtk_file%xml_reader%get_info(mesh_topology=topology)
    error = a_vtk_file%xml_reader%get_steps(times=times)
    print '(A,I0,A,*(F5.2))', repeat('  ', level(e))//'dataset '//trim(name(e))//' ('//trim(file(e))//', '// &
                              topology//'): ', size(times), ' steps, at', times
    error = a_vtk_file%finalize()
  endif
enddo

$ vtkhdf_motion
error: VTKFortran vtm_file: pier.vtkhdf is not a time series, the blocks before are (a file is either static or a time series)
the static block: error 1
harbour.vtkhdf, its blocks and datasets:
  block water
    dataset sea (Block_1, PolyData): 6 steps, at 0.00 0.25 0.50 0.75 1.00 1.25
  block bodies
    dataset buoy (Block_2, UnstructuredGrid): 6 steps, at 0.00 0.25 0.50 0.75 1.00 1.25

• A file is either static or a time series, as its first block: a static block in a time series is refused, as the output shows (VTK cannot read the mix). A part that does not move is written as a series too, the same mesh at every step.
• Every block has the same steps, their number and times; each one an UnstructuredGrid or a PolyData, of any number of parts per step. A time-series ImageData, RectilinearGrid or StructuredGrid block is refused: VTK 9.7.2 reads it without its arrays.
• ParaView 6.1.1 plays every step of every block. The rules and readers in detail: Multi-block in one VTKHDF file.

A time history (Table)

Time histories, probes and integrals are columns of values, without a mesh: a Table. Here the history of a damped oscillator, one row per time step, written at the end of the run as the columns of one Table:

! the history of x'' + 0.2 x' + x = 0, x(0) = 1, v(0) = 0 (semi-implicit Euler), one sample per time step
t(1) = 0 ; x(1) = 1 ; v(1) = 0
do i=2, n_samples
  t(i) = t(i-1) + dt
  v(i) = v(i-1) - dt * (0.2_R8P * v(i-1) + x(i-1))
  x(i) = x(i-1) + dt * v(i)
enddo
! one Table: a column per quantity, a row per sample
error = a_vtk_file%initialize(format='vtkhdf', filename='history.vtkhdf', mesh_topology='Table')
error = a_vtk_file%xml_writer%write_dataarray(location='row', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='time', x=t)
error = a_vtk_file%xml_writer%write_dataarray(data_name='displacement', x=x)
error = a_vtk_file%xml_writer%write_dataarray(data_name='velocity', x=v)
error = a_vtk_file%xml_writer%write_dataarray(location='row', action='close')
error = a_vtk_file%finalize()

A time series of Tables keeps the rows of each output apart: here three probes of a travelling wave, one row per probe at each output, the file valid after each step:

! a time series of Tables: at each output a row per probe, the same columns
probe_x = [1._R8P, 4._R8P, 7._R8P]
error = a_vtk_file%initialize(format='vtkhdf', filename='probes.vtkhdf', mesh_topology='Table')
do s=0, 3
  error = a_vtk_file%xml_writer%write_step(time=0.5_R8P * s)
  error = a_vtk_file%xml_writer%write_dataarray(location='row', action='open')
  error = a_vtk_file%xml_writer%write_dataarray(data_name='probe', x=[(i, i=1, n_probes)])
  error = a_vtk_file%xml_writer%write_dataarray(data_name='x', x=probe_x)
  error = a_vtk_file%xml_writer%write_dataarray(data_name='wave', x=sin(2 * pi * (probe_x - 0.5_R8P * s) / 8))
  error = a_vtk_file%xml_writer%write_dataarray(location='row', action='close')
enddo
error = a_vtk_file%finalize()

Both are read back with the reader of the library, the columns with location='row':

error = reader%initialize(filename='history.vtkhdf', action='read')
error = reader%xml_reader%read_dataarray(location='row', data_name='displacement', x=column)
print '(A,I0,A)', 'history.vtkhdf: ', size(column), ' rows; the displacement every 2 time units:'
print '(2X,*(F7.3))', column(1::8)
error = reader%finalize()
error = reader%initialize(filename='probes.vtkhdf', action='read')
error = reader%xml_reader%get_steps(times=times)
print '(A,I0,A)', 'probes.vtkhdf: ', size(times), ' steps; the wave at the probes:'
do s=1, size(times)
  error = reader%xml_reader%set_step(s)
  error = reader%xml_reader%read_dataarray(location='row', data_name='wave', x=column)
  print '(A,F4.2,A,*(F7.3))', '  t = ', times(s), ':', column
enddo
error = reader%finalize()

$ vtkhdf_table
history.vtkhdf: 41 rows; the displacement every 2 time units:
    1.000 -0.373 -0.399  0.530 -0.111 -0.273
probes.vtkhdf: 4 steps; the wave at the probes:
  t = 0.00:  0.707  0.000 -0.707
  t = 0.50:  0.383  0.383 -0.924
  t = 1.00:  0.000  0.707 -1.000
  t = 1.50: -0.383  0.924 -0.924

A Table is VTKHDF 2.8: VTK 9.7 reads it (vtkHDFReader gives a vtkTable), ParaView 6.1 does not (its VTK 9.6 does not know the type). See VTKHDF for the layout.

For ParaView, the same calls in an XML format write a .vtt (VTK's XML table, which ParaView opens as a spreadsheet and plots); a time series of Tables is then a .pvd of .vtt files:

! the history of x'' + 0.2 x' + x = 0, x(0) = 1, v(0) = 0 (semi-implicit Euler), one sample per time step
t(1) = 0 ; x(1) = 1 ; v(1) = 0
do i=2, n_samples
  t(i) = t(i-1) + dt
  v(i) = v(i-1) - dt * (0.2_R8P * v(i-1) + x(i-1))
  x(i) = x(i-1) + dt * v(i)
enddo
! one Table: a column per quantity, a row per sample (the same calls as for a VTKHDF Table)
error = a_vtk_file%initialize(format='binary', filename='history.vtt', mesh_topology='Table')
error = a_vtk_file%xml_writer%write_dataarray(location='row', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='time', x=t)
error = a_vtk_file%xml_writer%write_dataarray(data_name='displacement', x=x)
error = a_vtk_file%xml_writer%write_dataarray(data_name='velocity', x=v)
error = a_vtk_file%xml_writer%write_dataarray(location='row', action='close')
error = a_vtk_file%finalize()

! a time series of Tables: a .vtt per output, a row per probe, listed in a .pvd
probe_x = [1._R8P, 4._R8P, 7._R8P]
error = pvd%initialize(filename='probes.pvd')
do s=0, 3
  write(filename, '(A,I0,A)') 'probes_', s, '.vtt'
  error = a_vtk_file%initialize(format='ascii', filename=trim(filename), mesh_topology='Table')
  error = a_vtk_file%xml_writer%write_dataarray(location='row', action='open')
  error = a_vtk_file%xml_writer%write_dataarray(data_name='probe', x=[(i, i=1, n_probes)])
  error = a_vtk_file%xml_writer%write_dataarray(data_name='x', x=probe_x)
  error = a_vtk_file%xml_writer%write_dataarray(data_name='wave', x=sin(2 * pi * (probe_x - 0.5_R8P * s) / 8))
  error = a_vtk_file%xml_writer%write_dataarray(location='row', action='close')
  error = a_vtk_file%finalize()
  error = pvd%write_dataset(filename=trim(filename), timestep=0.5_R8P * s)
enddo
error = pvd%finalize()

error = reader%initialize(filename='history.vtt', action='read')
error = reader%xml_reader%get_dataarray_names(location='row', names=names)
print '(A,*(1X,A))', 'history.vtt columns:', (trim(names(i)), i=1, size(names))
error = reader%xml_reader%read_dataarray(location='row', data_name='displacement', x=column)
print '(A,I0,A)', 'history.vtt: ', size(column), ' rows; the displacement every 2 time units:'
print '(2X,*(F7.3))', column(1::8)
error = reader%finalize()

$ vtt_table
history.vtt columns: time displacement velocity
history.vtt: 41 rows; the displacement every 2 time units:
    1.000 -0.373 -0.399  0.530 -0.111 -0.273

The columns of a .vtt are kept in memory until write_dataarray(location='row', action='close'): its Piece tag, written first, holds their counts.

A block-structured AMR

A block-structured AMR in one file, an OverlappingAMR: level 0 covers the domain with one box, level 1 refines two boxes around a hot spot at half the spacing. Each level starts with write_level, each box is a piece of it (its point extent in the indices of the level), with its arrays:

error = a_vtk_file%initialize(format='vtkhdf', filename='amr.vtkhdf', mesh_topology='OverlappingAMR', &
                              origin=[0._R8P, 0._R8P, 0._R8P])
! level 0: the whole domain, spacing 1
error = a_vtk_file%xml_writer%write_level(level=0, spacing=[1._R8P, 1._R8P, 1._R8P])
error = box(nx1=0, nx2=16, ny1=0, ny2=8, h=1._R8P)
! level 1: spacing 0.5, two boxes around the hot spot (their point indices are those of level 1)
error = a_vtk_file%xml_writer%write_level(level=1, spacing=[0.5_R8P, 0.5_R8P, 0.5_R8P])
error = box(nx1=8, nx2=16, ny1=4, ny2=12, h=0.5_R8P)
error = box(nx1=16, nx2=24, ny1=4, ny2=12, h=0.5_R8P)
error = a_vtk_file%finalize()

The file is read back level by level, the boxes of a level as its pieces:

error = reader%initialize(filename='amr.vtkhdf', action='read')
error = reader%xml_reader%get_levels(n_levels)
do level=0, n_levels - 1
  error = reader%xml_reader%set_level(level)
  error = reader%xml_reader%get_info(npieces=n_boxes, spacing=spacing)
  print '(A,I0,A,F4.2,A,I0)', 'level ', level, ': spacing ', spacing(1), ', boxes ', n_boxes
  do b=1, n_boxes
    error = reader%xml_reader%read_piece(piece=b, nx1=nx1, nx2=nx2, ny1=ny1, ny2=ny2)
    error = reader%xml_reader%read_dataarray(location='node', data_name='temperature', x=temperature, piece=b)
    print '(A,I0,A,4(1X,I2),A,F6.3)', '  box ', b, ': points', nx1, nx2, ny1, ny2, ', hottest', maxval(temperature)
  enddo
enddo
error = reader%finalize()

$ vtkhdf_amr
level 0: spacing 1.00, boxes 1
  box 1: points  0 16  0  8, hottest 1.000
level 1: spacing 0.50, boxes 2
  box 1: points  8 16  4 12, hottest 0.368
  box 2: points 16 24  4 12, hottest 1.000

VTK 9.7 reads it (vtkHDFReader gives a vtkOverlappingAMR); ParaView 6.1 does not: its VTKHDF reader passes the output to a filter that refuses an OverlappingAMR. See VTKHDF for the layout.

A tree-based AMR

A tree-based AMR in one file, a HyperTreeGrid: a root grid of 4x2 cells, each the root of a quadtree refined down to depth 3 around a circle. Each tree is its refinement, breadth first (the deepest depth left out: its cells are all leaves); the depth of each cell is its cell data:

error = a_vtk_file%initialize(format='vtkhdf', filename='quadtree.vtkhdf', mesh_topology='HyperTreeGrid', &
                              nx1=0, nx2=nx, ny1=0, ny2=ny, nz1=0, nz2=0, branch_factor=2)
error = a_vtk_file%xml_writer%write_piece(nx1=0, nx2=nx, ny1=0, ny2=ny, nz1=0, nz2=0)
error = a_vtk_file%xml_writer%write_geo(x=[(real(i, R8P), i=0, nx)], y=[(real(j, R8P), j=0, ny)], z=[0._R8P])
allocate(depths(0))
do j=0, ny - 1 ; do i=0, nx - 1 ! the tree of root cell (i, j): its id i + nx j
  error = a_vtk_file%xml_writer%write_tree(tree_id=i + nx * j, descriptor=tree(real(i, R8P), real(j, R8P)))
enddo ; enddo
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='depth', x=depths)
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

Read back, each tree comes as write_tree took it:

error = reader%initialize(filename='quadtree.vtkhdf', action='read')
error = reader%xml_reader%get_trees(tree_ids=ids, branch_factor=bf)
print '(A,I0,A,I0)', 'read back: ', size(ids), ' trees, branch factor ', bf
do t=1, size(ids) ! each tree as write_tree took it
  error = reader%xml_reader%read_tree(tree_id=ids(t), descriptor=descriptor)
  print '(A,I0,A,I0,A,I0,A)', '  tree ', ids(t), ': ', size(descriptor), ' descriptor bits, ', count(descriptor), ' refined'
enddo
error = reader%xml_reader%read_dataarray(location='cell', data_name='depth', x=depths_read)
print '(A,L1)', '  the cell data as written: ', all(depths_read == depths)
error = reader%finalize()

$ vtkhdf_htg
quadtree.vtkhdf: 8 trees, 184 cells
  cells of each depth: 8 16 48 112
read back: 8 trees, branch factor 2
  tree 0: 0 descriptor bits, 0 refined
  tree 1: 17 descriptor bits, 11 refined
  tree 2: 17 descriptor bits, 11 refined
  tree 3: 0 descriptor bits, 0 refined
  tree 4: 0 descriptor bits, 0 refined
  tree 5: 17 descriptor bits, 11 refined
  tree 6: 17 descriptor bits, 11 refined
  tree 7: 0 descriptor bits, 0 refined
  the cell data as written: T

The program builds each tree level by level (the function tree of vtkhdf_htg.f90). VTK 9.7 and ParaView 6.1 read it. See VTKHDF for the layout.

More than 2^31 points or cells

program ids_64bit
!< Write the counts and the connectivity in 64 bits, for meshes with more than 2^31 points or cells.
use penf, only : I1P, I8P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file) :: a_vtk_file
integer        :: error

! the kind of the arguments selects the version: I8P counts and ids are written as Int64
error = a_vtk_file%initialize(format='ascii', filename='tet64.vtu', mesh_topology='UnstructuredGrid')
error = a_vtk_file%xml_writer%write_piece(np=4_I8P, nc=1_I8P)
error = a_vtk_file%xml_writer%write_geo(np=4_I8P, nc=1_I8P, x=[0._R8P, 1._R8P, 0._R8P, 0._R8P], &
                                        y=[0._R8P, 0._R8P, 1._R8P, 0._R8P], z=[0._R8P, 0._R8P, 0._R8P, 1._R8P])
error = a_vtk_file%xml_writer%write_connectivity(nc=1_I8P, connectivity=[0_I8P, 1_I8P, 2_I8P, 3_I8P], offset=[4_I8P], &
                                                 cell_type=[10_I1P])
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
endprogram ids_64bit

$ grep -o 'type="Int64"[^>]*Name="[a-z]*"' tet64.vtu
type="Int64" NumberOfComponents="1" Name="connectivity"
type="Int64" NumberOfComponents="1" Name="offsets"

The kind of the counts and ids selects the version: I4P calls are unchanged.

Several pieces in one file

program pieces_in_file
!< Write several pieces in one file: e.g. the blocks of a solver, written by one process.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file) :: a_vtk_file
integer(I4P)   :: b, error

! the whole extent, then each piece with its own extent inside it; adjacent pieces share their boundary
error = a_vtk_file%initialize(format='raw', filename='blocks.vtr', mesh_topology='RectilinearGrid', &
                              nx1=0, nx2=6, ny1=0, ny2=1, nz1=0, nz2=1)
do b=1, 3
  error = a_vtk_file%xml_writer%write_piece(nx1=2*(b-1), nx2=2*b, ny1=0, ny2=1, nz1=0, nz2=1)
  error = a_vtk_file%xml_writer%write_geo(x=[real(2*b-2, R8P), real(2*b-1, R8P), real(2*b, R8P)], y=[0._R8P, 1._R8P], &
                                          z=[0._R8P, 1._R8P])
  error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
  error = a_vtk_file%xml_writer%write_dataarray(data_name='block', x=[b, b])
  error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
  error = a_vtk_file%xml_writer%write_piece()
enddo
error = a_vtk_file%finalize()
endprogram pieces_in_file

$ inspect blocks.vtr
blocks.vtr: RectilinearGrid, pieces: 3, compressor: none
  piece 1: 12 points, 2 cells
     cell block              Int32 1 comp., range  1.0000E+00  1.0000E+00
  piece 2: 12 points, 2 cells
     cell block              Int32 1 comp., range  2.0000E+00  2.0000E+00
  piece 3: 12 points, 2 cells
     cell block              Int32 1 comp., range  3.0000E+00  3.0000E+00

Readers merge the pieces into one dataset. Every piece must hold the same arrays.

Parallel and composite files

A partitioned structured grid (.pvts)

program pvts
!< Write a structured grid split in two pieces, as two processes would, and the .pvts header that lists them.
use penf, only : I4P, R8P
use vtk_fortran, only : pvtk_file, vtk_file
implicit none
type(pvtk_file)               :: header
character(len=:), allocatable :: message
integer(I4P)                  :: error

! each process: a complete .vts with the extent of its piece (adjacent pieces share the plane i=3)
call write_piece('part_1.vts', 1, 3)
call write_piece('part_2.vts', 3, 5)
! one process: the header, with the layout of the data and the extent of each piece
error = header%initialize(filename='grid.pvts', mesh_topology='PStructuredGrid', mesh_kind='Float64', &
                          nx1=1, nx2=5, ny1=1, ny2=2, nz1=1, nz2=2)
error = header%xml_writer%write_dataarray(location='node', action='open')
error = header%xml_writer%write_parallel_dataarray(data_name='x', data_type='Float64', number_of_components=1)
error = header%xml_writer%write_dataarray(location='node', action='close')
error = header%xml_writer%write_parallel_geo(source='part_1.vts', nx1=1, nx2=3, ny1=1, ny2=2, nz1=1, nz2=2)
error = header%xml_writer%write_parallel_geo(source='part_2.vts', nx1=3, nx2=5, ny1=1, ny2=2, nz1=1, nz2=2)
error = header%finalize()
! check the pieces against the header
error = header%initialize(filename='grid.pvts', action='read')
error = header%xml_reader%check_pieces(message=message)
print '(A,I0)', 'grid.pvts: check_pieces error ', error
error = header%finalize()
contains
  subroutine write_piece(filename, i1, i2)
  !< Write the piece of the points i1...i2 along x.
  character(*), intent(in) :: filename
  integer(I4P), intent(in) :: i1, i2
  real(R8P)                :: x(i1:i2,2,2), y(i1:i2,2,2), z(i1:i2,2,2)
  type(vtk_file)           :: a_vtk_file
  integer(I4P)             :: i, error

  do i=i1, i2
    x(i,:,:) = i ; y(i,1,:) = 0 ; y(i,2,:) = 1 ; z(i,:,1) = 0 ; z(i,:,2) = 1
  enddo
  error = a_vtk_file%initialize(format='raw', filename=filename, mesh_topology='StructuredGrid', &
                                nx1=1, nx2=5, ny1=1, ny2=2, nz1=1, nz2=2)
  error = a_vtk_file%xml_writer%write_piece(nx1=i1, nx2=i2, ny1=1, ny2=2, nz1=1, nz2=2)
  error = a_vtk_file%xml_writer%write_geo(n=size(x), x=x, y=y, z=z)
  error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
  error = a_vtk_file%xml_writer%write_dataarray(data_name='x', x=x, one_component=.true.)
  error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
  error = a_vtk_file%xml_writer%write_piece()
  error = a_vtk_file%finalize()
  endsubroutine write_piece
endprogram pvts

$ pvts
grid.pvts: check_pieces error 0

The same for .pvtr (PRectilinearGrid) and .pvti (PImageData, with origin and spacing); unstructured pieces (.pvtu, .pvtp) have no extent, see chapter 6.

Write from MPI ranks

Each rank writes its own piece, a complete file named after the rank; rank 0 waits for all of them, writes the header that lists them and checks it. No data move between the ranks: every file is written by one process.

program mpi_pieces
!< Write a partitioned mesh from MPI ranks: each rank writes its own piece, rank 0 the parallel header, then checks it.
use mpi_f08
use penf, only : I1P, I4P, R8P
use vtk_fortran, only : pvtk_file, vtk_file
implicit none
integer, parameter            :: nc=6               ! cells of each piece: a strip of tetrahedra
integer                       :: rank, nranks, r
character(len=32)             :: filename
character(len=:), allocatable :: message
type(pvtk_file)               :: header
integer(I4P)                  :: error

call MPI_Init()
call MPI_Comm_rank(MPI_COMM_WORLD, rank)
call MPI_Comm_size(MPI_COMM_WORLD, nranks)

! every rank: its own piece, in a file named after the rank, with point ids local to the piece
write(filename, '(A,I0,A)') 'mpi_piece_', rank, '.vtu'
error = write_piece(trim(filename), rank)

call MPI_Barrier(MPI_COMM_WORLD)    ! the header is checked against the pieces: wait until all are written
if (rank == 0) then
  error = header%initialize(filename='mpi.pvtu', mesh_topology='PUnstructuredGrid', mesh_kind='Float64')
  error = header%xml_writer%write_dataarray(location='cell', action='open')
  error = header%xml_writer%write_parallel_dataarray(data_name='rank', data_type='Int32', number_of_components=1)
  error = header%xml_writer%write_dataarray(location='cell', action='close')
  do r=0, nranks - 1
    write(filename, '(A,I0,A)') 'mpi_piece_', r, '.vtu'
    error = header%xml_writer%write_parallel_geo(source=trim(filename))
  enddo
  error = header%finalize()
  ! read the header back and check every piece against it
  error = header%initialize(filename='mpi.pvtu', action='read')
  error = header%xml_reader%check_pieces(message=message)
  print '(A,I0,A,I0)', 'mpi.pvtu: ', nranks, ' pieces, check_pieces error ', error
  error = header%finalize()
endif
call MPI_Finalize()
contains
  function write_piece(filename, rank) result(error)
  !< Write the piece of a rank: a strip of nc tetrahedra, shifted along x by the rank, with the cell data "rank".
  character(*), intent(in) :: filename
  integer,      intent(in) :: rank
  integer(I4P)             :: error
  type(vtk_file)           :: a_vtk_file
  real(R8P)                :: x(nc+3), y(nc+3), z(nc+3)
  integer(I4P)             :: connect(4*nc), c, p

  do p=1, nc + 3
    x(p) = rank*nc + p - 1 ; y(p) = mod(p, 2) ; z(p) = mod(p, 3)
  enddo
  connect = [((c - 1 + p, p=0, 3), c=1, nc)]
  error = a_vtk_file%initialize(format='raw', filename=filename, mesh_topology='UnstructuredGrid')
  error = a_vtk_file%xml_writer%write_piece(np=nc+3, nc=nc)
  error = a_vtk_file%xml_writer%write_geo(np=nc+3, nc=nc, x=x, y=y, z=z)
  error = a_vtk_file%xml_writer%write_connectivity(nc=nc, connectivity=connect, offset=[(4*c, c=1, nc)], &
                                                   cell_type=[(10_I1P, c=1, nc)])
  error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
  error = a_vtk_file%xml_writer%write_dataarray(data_name='rank', x=[(int(rank, I4P), c=1, nc)])
  error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
  error = a_vtk_file%xml_writer%write_piece()
  error = a_vtk_file%finalize()
  endfunction write_piece
endprogram mpi_pieces

$ mpirun -np 4 mpi_pieces
mpi.pvtu: 4 pieces, check_pieces error 0

Built with the MPI compiler wrapper (mpifort) and run with mpirun -np 4. When only rank 0 can write to the file system, the other ranks write their pieces into memory and send them to it, see Write a file into memory.

Write one .vtkhdf from MPI ranks

With a parallel HDF5 the ranks write one file together: comm (the Fortran handle of a communicator) given to initialize, every rank makes the same calls, each one with its own piece, which becomes a partition of the file (an UnstructuredGrid, a PolyData), or its own region of the grid (ImageData, RectilinearGrid, StructuredGrid, below).

! every rank makes the same calls (they are collective), each one with its own piece
error = a_vtk_file%initialize(format='vtkhdf', filename='ranks.vtkhdf', mesh_topology='UnstructuredGrid', &
                              comm=MPI_COMM_WORLD%MPI_VAL)
call block_mesh(rank, x, y, z, connect, offset, cell_type) ! the hexahedra of this rank, point ids from 0
error = a_vtk_file%xml_writer%write_piece(np=np, nc=nc)
error = a_vtk_file%xml_writer%write_geo(np=np, nc=nc, x=x, y=y, z=z)
error = a_vtk_file%xml_writer%write_connectivity(nc=nc, connectivity=connect, offset=offset, cell_type=cell_type)
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='rank', x=[(rank, c=1, nc)])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

$ mpirun -np 4 vtkhdf_mpi
ranks.vtkhdf, one partition per rank:
  NumberOfPoints 18 18 18 18
  NumberOfCells  4 4 4 4
  CellData/rank  0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3

• Every call is collective: the ranks create and extend the datasets together, then each one writes its rows after those of the ranks before. Every rank calls every procedure, with the same arrays in the same order. Wrong arguments on some ranks in a call they all make (a size, a negative count, face without faceoffset) are refused on every rank, and so are ranks writing different arrays in the same call; but a rank leaving a call out (one stops calling, or, in a step reusing the mesh, writes an array the others do not) cannot be told apart from a slow one: the others wait for it.
• An UnstructuredGrid or a PolyData, one piece per rank; a rank can have an empty piece (np=0, nc=0). A PolyData piece without cells still calls write_polydata_cells(), with no blocks.
• Polyhedra (on any ranks: the others write zero entries for their pieces) and compressor='zlib' (an HDF5 1.10.2 or later with the deflate filter: a deflated dataset needs collective writes, which the ranks then make) are written as in serial. Refused with comm: compressor='zlib' for a grid, the XML formats, action='read'. A refusal comes on every rank, also when only one rank's values are wrong (their size, say): the ranks agree on the error before the collective calls, none is left waiting.

A time series from the ranks: every rank calls write_step with the same time (checked: different times are an error on every rank), then writes its piece of the step; a step without pieces reuses the mesh of the step before, each rank writing the arrays of its own part. Here the blocks of the four ranks, the mesh written at the first step only:

! every rank makes the same calls (they are collective): the same time, its own piece or point data
error = a_vtk_file%initialize(format='vtkhdf', filename='series.vtkhdf', mesh_topology='UnstructuredGrid', &
                              comm=MPI_COMM_WORLD%MPI_VAL)
do s=0, 5
  error = a_vtk_file%xml_writer%write_step(time=0.25_R8P * s)
  if (s == 0) then ! the first step: the mesh of each rank
    error = a_vtk_file%xml_writer%write_piece(np=np, nc=nc)
    error = a_vtk_file%xml_writer%write_geo(np=np, nc=nc, x=x, y=y, z=z)
    error = a_vtk_file%xml_writer%write_connectivity(nc=nc, connectivity=connect, offset=offset, cell_type=cell_type)
    error = point_data(time=0._R8P)
    error = a_vtk_file%xml_writer%write_piece()
  else             ! the next steps: no piece, the point data of each rank's block (the mesh of the step before)
    error = point_data(time=0.25_R8P * s)
  endif
enddo
error = a_vtk_file%finalize()

$ mpirun -np 4 vtkhdf_mpi_series
series.vtkhdf, the steps:
  NumberOfParts            4   4   4   4   4   4
  PointOffsets             0   0   0   0   0   0
  PointDataOffsets/wave    0  72 144 216 288 360

• The Steps entries are those of all the ranks: a part per rank in each step (NumberOfParts), offsets of the totals of the ranks before, written once. A rank's arrays in a step without pieces hold the points (cells) of its own part.
• initialize(format='vtkhdf', filename=..., action='append', comm=...) continues the series from the ranks (a parallel restart); a step reusing the mesh after it needs the same number of ranks as the parts of the last step.
• A time series of a grid works the same way, a region per rank in each step: see below.

A grid (ImageData, RectilinearGrid, StructuredGrid) is written by regions: initialize gets the whole extent, write_piece the region of the rank, and its arrays (and coordinates or points) are those of the region. Here an ImageData of 8x4 cells, one region of 4x2 cells per rank:

! the region of this rank: 4x2 cells of the whole extent 0:8 x 0:4 (2D)
nx1 = 4 * mod(rank, 2) ; nx2 = nx1 + 4 ; ny1 = 2 * (rank / 2) ; ny2 = ny1 + 2
wave = [((sin(0.4_R8P * i) * cos(0.6_R8P * j), i=nx1, nx2), j=ny1, ny2)] ! its points, global indices
! every rank makes the same calls (they are collective): the whole extent to initialize, its region to write_piece
error = a_vtk_file%initialize(format='vtkhdf', filename='grid.vtkhdf', mesh_topology='ImageData', &
                              nx1=0, nx2=8, ny1=0, ny2=4, nz1=0, nz2=0,                          &
                              origin=[0._R8P, 0._R8P, 0._R8P], spacing=[1._R8P, 1._R8P, 1._R8P], &
                              comm=MPI_COMM_WORLD%MPI_VAL)
error = a_vtk_file%xml_writer%write_piece(nx1=nx1, nx2=nx2, ny1=ny1, ny2=ny2, nz1=0, nz2=0)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='wave', x=wave)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='close')
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='rank', x=[(rank, i=1, 8)])
error = a_vtk_file%xml_writer%write_dataarray(location='cell', action='close')
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()

$ mpirun -np 4 vtkhdf_mpi_grid
grid.vtkhdf: PointData/wave, Fortran shape 9 5
CellData/rank, the 8x4 cells (the row y=0 last):
   2 2 2 2 3 3 3 3
   2 2 2 2 3 3 3 3
   0 0 0 0 1 1 1 1
   0 0 0 0 1 1 1 1

• The regions of all the ranks are checked when they start, on every rank: inside the whole extent, not overlapping in cells, covering it. A rank with nothing to write gives an empty region (an upper bound before its lower one) and empty arrays.
• Neighbouring regions share the points of their boundary: each rank gives all the points of its region, and the boundary is written by the region that starts there. The coordinates of a RectilinearGrid shared by regions must be equal; the lowest rank writes them.
• The file is the one the whole grid written by one process would make: VTK and ParaView read it as any ImageData.
• A time series of a grid: every rank calls write_step with the same time, then write_piece with its region of the step (checked at every step: it may change, the whole extent does not) and the arrays of the region. Each dataset grows by one time slice per step, each rank writing its block of it; a step without write_geo reuses the geometry of the step before, an array not written its values (the ranks make the same calls). action='append' with comm continues the series (a parallel restart): the first step after it writes its piece (the regions are not in the file; a step without pieces there is refused). The file is, dataset by dataset, the one the whole grid written by one process would make (src/tests/vtk_fortran_write_vtkhdf_ranks.F90 checks it on 1 to 4 ranks, VTK 9.7.2 each step against its XML file).
• It needs the library built with VTKFORTRAN_USE_MPI (CMake -DVTKFORTRAN_USE_MPI=ON, FoBiS mode static-gnu-hdf5-mpi) and a parallel HDF5; initialize returns an error otherwise (Installation). Built with mpifort, run with mpirun -np 4.

Check a parallel header against its pieces

program check_pieces
!< Find the mismatch between a parallel header and its pieces before ParaView does.
use penf, only : I1P, I4P, R4P, R8P
use vtk_fortran, only : pvtk_file, vtk_file
implicit none
type(pvtk_file)               :: header
type(vtk_file)                :: piece
character(len=:), allocatable :: message
integer(I4P)                  :: error

! a piece with a Float32 pressure
error = piece%initialize(format='raw', filename='p_1.vtu', mesh_topology='UnstructuredGrid')
error = piece%xml_writer%write_piece(np=4, nc=1)
error = piece%xml_writer%write_geo(np=4, nc=1, x=[0._R8P, 1._R8P, 0._R8P, 0._R8P], y=[0._R8P, 0._R8P, 1._R8P, 0._R8P], &
                                   z=[0._R8P, 0._R8P, 0._R8P, 1._R8P])
error = piece%xml_writer%write_connectivity(nc=1, connectivity=[0, 1, 2, 3], offset=[4], cell_type=[10_I1P])
error = piece%xml_writer%write_dataarray(location='node', action='open')
error = piece%xml_writer%write_dataarray(data_name='pressure', x=[1._R4P, 2._R4P, 3._R4P, 4._R4P])
error = piece%xml_writer%write_dataarray(location='node', action='close')
error = piece%xml_writer%write_piece()
error = piece%finalize()
! a header declaring it Float64
error = header%initialize(filename='p.pvtu', mesh_topology='PUnstructuredGrid', mesh_kind='Float64')
error = header%xml_writer%write_dataarray(location='node', action='open')
error = header%xml_writer%write_parallel_dataarray(data_name='pressure', data_type='Float64', number_of_components=1)
error = header%xml_writer%write_dataarray(location='node', action='close')
error = header%xml_writer%write_parallel_geo(source='p_1.vtu')
error = header%finalize()

error = header%initialize(filename='p.pvtu', action='read')
error = header%xml_reader%check_pieces(message=message)
print '(A,I0)', 'check_pieces error ', error
print '(A)', message
error = header%finalize()
endprogram check_pieces

$ check_pieces
check_pieces error 7
piece 1 (p_1.vtu): node array "pressure" is Float32 with 1 components, declared Float64 with 1

Continue a time series after a restart

program pvd_append
!< Continue a time series after a restart, then read the collection back.
use penf, only : I4P, R8P
use vtk_fortran, only : pvd_file
implicit none
type(pvd_file)                :: series
real(R8P),        allocatable :: timestep(:)
character(len=:), allocatable :: file(:)
integer(I4P)                  :: d, error

! the first job
error = series%initialize(filename='run.pvd')
error = series%write_dataset(filename='run_0.vtu', timestep=0._R8P)
error = series%write_dataset(filename='run_1.vtu', timestep=0.5_R8P)
error = series%finalize()
! the restarted job: the datasets already listed are kept
error = series%initialize(filename='run.pvd', action='append')
error = series%write_dataset(filename='run_2.vtu', timestep=1._R8P)
error = series%finalize()

error = series%initialize(filename='run.pvd', action='read')
error = series%get_datasets(timestep=timestep, file=file)
error = series%finalize()
do d=1, size(file)
  print '(A,F4.2,2A)', 'time ', timestep(d), ': ', file(d)
enddo
endprogram pvd_append

$ pvd_append
time 0.00: run_0.vtu
time 0.50: run_1.vtu
time 1.00: run_2.vtu

The collection is valid after every write_dataset: a job killed before finalize leaves a readable series.

Assemblies of datasets (.vtm)

Blocks nested to any depth, read back with get_entries: see chapter 7.

Write a file into memory

program volatile
!< Write a file into memory, as a process without access to the file system would; another one saves it.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file, write_xml_volatile
implicit none
type(vtk_file)                :: a_vtk_file
character(len=:), allocatable :: xml
integer(I4P)                  :: error

error = a_vtk_file%initialize(format='binary', filename='part.vtr', mesh_topology='RectilinearGrid', &
                              nx1=1, nx2=2, ny1=1, ny2=2, nz1=1, nz2=2, is_volatile=.true.)
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=2, ny1=1, ny2=2, nz1=1, nz2=2)
error = a_vtk_file%xml_writer%write_geo(x=[0._R8P, 1._R8P], y=[0._R8P, 1._R8P], z=[0._R8P, 1._R8P])
error = a_vtk_file%xml_writer%write_piece()
error = a_vtk_file%finalize()
call a_vtk_file%get_xml_volatile(xml) ! the whole file, as a string: send it to the process that writes
call a_vtk_file%free
print '(A,I0,A)', 'the file is in memory: ', len(xml), ' characters'
! ... on the process that accesses the file system
error = write_xml_volatile(xml_volatile=xml, filename='part.vtr')
print '(A,I0)', 'part.vtr written, error ', error
endprogram volatile

$ volatile
the file is in memory: 722 characters
part.vtr written, error 0

For the processes that cannot access the file system: the binary and ascii formats support volatile files; the appended ones (raw, raw-zlib, binary-appended) write their data to the file directly, and initialize refuses them.

Read files

What a file holds

inspect, the program of the quick start, lists the pieces and the arrays of any file, with their range.

Read one array

program read_array
!< Read one array of a file (here a file written by VTK): ask its type and size, then read it into a kind that holds it.
use penf, only : I4P, I8P, R4P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file)                :: a_vtk_file
character(len=:), allocatable :: data_type
real(R8P),        allocatable :: p(:)
integer(I4P),     allocatable :: wrong(:)
integer(I8P)                  :: n_tuples
integer(I4P)                  :: n_components, error

error = a_vtk_file%initialize(filename='vtk_tetra.vtu', action='read')
error = a_vtk_file%xml_reader%get_dataarray_info(location='node', data_name='p', data_type=data_type, &
                                                 n_components=n_components, n_tuples=n_tuples)
print '(A,A,A,I0,A,I0,A)', 'p: ', data_type, ', ', n_components, ' component, ', n_tuples, ' tuples'
! Float32 reads into R4P or R8P
error = a_vtk_file%xml_reader%read_dataarray(location='node', data_name='p', x=p)
print '(A,4F5.1)', 'p:', p
! not into an integer: error 5
error = a_vtk_file%xml_reader%read_dataarray(location='node', data_name='p', x=wrong)
print '(A,I0)', 'p into I4P: error ', error
error = a_vtk_file%finalize()
endprogram read_array

$ read_array
p: Float32, 1 component, 4 tuples
p:  1.5  2.5  3.5  4.5
p into I4P: error 5

The output kind must hold every value of the type: an error 5 otherwise. Unsigned types read into the signed kind of the same width, with the same bits.

Read the mesh

program read_mesh
!< Read the points and the cells of an unstructured grid (here written by VTK, with Int64 ids and zlib compression).
use penf, only : I1P, I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file)            :: a_vtk_file
real(R8P),    allocatable :: x(:), y(:), z(:)
integer(I4P), allocatable :: connectivity(:), offset(:)
integer(I1P), allocatable :: cell_type(:)
integer(I4P)              :: i, error

error = a_vtk_file%initialize(filename='vtk_tetra.vtu', action='read')
error = a_vtk_file%xml_reader%read_geo(x=x, y=y, z=z)
! the ids are Int64 in the file: they read into I4P, since they fit
error = a_vtk_file%xml_reader%read_connectivity(connectivity=connectivity, offset=offset, cell_type=cell_type)
error = a_vtk_file%finalize()
do i=1, size(x)
  print '(A,I0,A,3F5.1)', 'point ', i - 1, ':', x(i), y(i), z(i)
enddo
print '(A,I0)', 'cell type ', cell_type(1)
print '(A,*(1X,I0))', 'its points', connectivity
endprogram read_mesh

$ read_mesh
point 0:  0.0  0.0  0.0
point 1:  1.0  0.0  0.0
point 2:  0.0  1.0  0.0
point 3:  0.0  0.0  1.0
cell type 10
its points 0 1 2 3

vtk_tetra.vtu was written by VTK (zlib compressed, base64 appended, UInt64 headers): the reader takes the format from the file.

Handle the errors

program errors
!< Every procedure returns an error status: 0 on success; the readers say what went wrong.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
type(vtk_file)         :: a_vtk_file
real(R8P), allocatable :: x(:)
integer(I4P)           :: error

error = a_vtk_file%initialize(filename='no_such_file.vtu', action='read')
print '(A,I0)', 'a missing file:        ', error  ! 1
error = a_vtk_file%initialize(filename='vtk_tetra.vtu', action='read')
print '(A,I0)', 'a good file:           ', error  ! 0
error = a_vtk_file%xml_reader%read_dataarray(location='node', data_name='no_such_array', x=x)
print '(A,I0)', 'a missing array:       ', error  ! 4
error = a_vtk_file%xml_reader%read_dataarray(location='node', data_name='p', x=x, piece=2)
print '(A,I0)', 'a missing piece:       ', error  ! 4
error = a_vtk_file%finalize()
error = a_vtk_file%initialize(format='rawest', filename='x.vtu', mesh_topology='UnstructuredGrid')
print '(A,I0)', 'an unknown format:     ', error  ! not 0
endprogram errors

$ errors
a missing file:        1
a good file:           0
a missing array:       4
a missing piece:       4
an unknown format:     1

The codes of the readers: 1 the file cannot be read, 2 not a VTK XML file, 3 an unsupported feature, 4 not found, 5 a kind that cannot hold the values, 6 data that do not decode, 7 pieces that do not match their header.