Skip to content

1. A first file ​

heat starts from the initial temperature: two hot blobs in a cold cube, sampled at 24 points along each side. The points are aligned along the axes, so the cube is a rectilinear grid: the coordinates along each axis are enough to define all the points.

f90
x = [(real(i - 1, R8P)/(n - 1), i=1, n)]
do k=1, n ; do j=1, n ; do i=1, n
  ! two hot blobs in a cold cube: the walls are kept at 0
  t(i,j,k) = blob(x(i), x(j), x(k), [0.35_R8P, 0.4_R8P, 0.5_R8P]) + 0.6_R8P*blob(x(i), x(j), x(k), [0.7_R8P, 0.65_R8P, 0.45_R8P])
enddo ; enddo ; enddo
t(1,:,:) = 0 ; t(n,:,:) = 0 ; t(:,1,:) = 0 ; t(:,n,:) = 0 ; t(:,:,1) = 0 ; t(:,:,n) = 0

Writing it takes one call for each part of the file:

f90
error = a_vtk_file%initialize(format='ascii', filename='heat.vtr', mesh_topology='RectilinearGrid', &
                              nx1=1, nx2=n, ny1=1, ny2=n, nz1=1, nz2=n)
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=n, ny1=1, ny2=n, nz1=1, nz2=n)
error = a_vtk_file%xml_writer%write_geo(x=x, y=x, z=x)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='temperature', x=t, 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()
  • initialize opens the file: the format (ascii here, human readable), the file name, the type of the dataset (mesh_topology) and its whole extent, the range of the point indexes along each axis.
  • A file holds one or more pieces: write_piece with the extent opens one, write_piece() without arguments closes it.
  • write_geo writes the geometry; for a rectilinear grid, the coordinates along each axis.
  • The point data are written between an open and a close of the node location; write_dataarray accepts arrays of rank 1 to 4: one_component=.true. says that the rank-3 t is one scalar per point, not a vector.
  • finalize closes the file. Every procedure returns an error status: 0 means success.

The whole program:

heat_1.f90
f90
program heat
!< Tutorial, chapter 1: the initial temperature of the cube, written as a rectilinear grid in ASCII.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file
implicit none
integer(I4P), parameter :: n=24          ! points along each side of the unit cube
real(R8P)               :: x(n)          ! coordinates of the points along each axis
real(R8P)               :: t(n,n,n)      ! temperature at the points
type(vtk_file)          :: a_vtk_file
integer(I4P)            :: i, j, k, error

x = [(real(i - 1, R8P)/(n - 1), i=1, n)]
do k=1, n ; do j=1, n ; do i=1, n
  ! two hot blobs in a cold cube: the walls are kept at 0
  t(i,j,k) = blob(x(i), x(j), x(k), [0.35_R8P, 0.4_R8P, 0.5_R8P]) + 0.6_R8P*blob(x(i), x(j), x(k), [0.7_R8P, 0.65_R8P, 0.45_R8P])
enddo ; enddo ; enddo
t(1,:,:) = 0 ; t(n,:,:) = 0 ; t(:,1,:) = 0 ; t(:,n,:) = 0 ; t(:,:,1) = 0 ; t(:,:,n) = 0

error = a_vtk_file%initialize(format='ascii', filename='heat.vtr', mesh_topology='RectilinearGrid', &
                              nx1=1, nx2=n, ny1=1, ny2=n, nz1=1, nz2=n)
error = a_vtk_file%xml_writer%write_piece(nx1=1, nx2=n, ny1=1, ny2=n, nz1=1, nz2=n)
error = a_vtk_file%xml_writer%write_geo(x=x, y=x, z=x)
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='temperature', x=t, 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,A,F6.3)', 'heat.vtr: ', n**3, ' points, maximum temperature ', maxval(t)
contains
  pure function blob(x, y, z, centre) result(t)
  !< A hot blob: a Gaussian bump of temperature 1 at its centre.
  real(R8P), intent(in) :: x, y, z, centre(3)
  real(R8P)             :: t

  t = exp(-((x - centre(1))**2 + (y - centre(2))**2 + (z - centre(3))**2)/0.04_R8P)
  endfunction blob
endprogram heat

Running it ​

$ heat
heat.vtr: 13824 points, maximum temperature  0.991

The file is plain XML; its head:

$ head -n 6 heat.vtr
<?xml version="1.0"?>
<VTKFile type="RectilinearGrid" version="1.0" byte_order="LittleEndian" header_type="UInt32">
  <RectilinearGrid WholeExtent="+1 +24 +1 +24 +1 +24">
    <Piece Extent="+1 +24 +1 +24 +1 +24">
      <Coordinates>
        <DataArray type="Float64" NumberOfComponents="1" Name="X" format="ascii">

Opened in ParaView (here with a few isosurfaces of the temperature, cut in half):

isosurfaces of the temperature of two hot blobs in a cube

What you learned

A file is initialize, then pieces with their geometry and data, then finalize. Reference: Rectilinear Grid.

Next: 2. Formats.