Appearance
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)
No coordinates at all: the points are defined by the extent, an origin and a spacing.
f90
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.
f90
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.
f90
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
integer(I4P), parameter :: face(31)=[6, 4,0,1,2,3, 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.
f90
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.
f90
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+01scalars, vectors and tensors choose the arrays readers use by default. inspect is the program of the quick start.
Time, cycle and other global data
f90
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 CO2Field 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
f90
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 1write_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.
Choose the format
Compress, and go beyond 2 GiB
f90
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 zlibcompressor='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.
More than 2^31 points or cells
f90
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
f90
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+00Readers merge the pieces into one dataset. Every piece must hold the same arrays.
Parallel and composite files
A partitioned structured grid (.pvts)
f90
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 0The same for .pvtr (PRectilinearGrid) and .pvti (PImageData, with origin and spacing); unstructured pieces (.pvtu, .pvtp) have no extent, see chapter 6.
Check a parallel header against its pieces
f90
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 1Continue a time series after a restart
f90
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.vtuThe 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
f90
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 0For 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
f90
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 5The 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
f90
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 3vtk_tetra.vtu was written by VTK (zlib compressed, base64 appended, UInt64 headers): the reader takes the format from the file.
Handle the errors
f90
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: 1The 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.