7. An assembly

A simulation often has more than one dataset: the domain, the probes that record the temperature at a few points, a boundary, a body. A multi-block .vtm file gathers them in blocks, nested as you like, and ParaView opens it as one object whose parts can be shown or hidden.

The probes are a polydata of 8 points, one vertex each, with the temperature read from the pieces of chapter 6:

! the probes: a polydata of points, one vertex each
error = a_vtk_file%initialize(format='ascii', filename='probes.vtp', mesh_topology='PolyData')
error = a_vtk_file%xml_writer%write_piece(np=probes, nverts=probes, nlines=0, nstrips=0, npolys=0)
error = a_vtk_file%xml_writer%write_geo(np=probes, nc=probes, x=px, y=py, z=pz)
error = a_vtk_file%xml_writer%write_polydata_cells(verts_connectivity=[(e - 1, e=1, probes)], verts_offset=[(e, e=1, probes)])
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='temperature', x=tp)
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 assembly: a solver block holding the domain (the 4 pieces) and the sensors:

error = assembly%initialize(filename='heat.vtm')
error = assembly%write_block(action='open', name='solver')
error = assembly%write_block(filenames=['heat_1.vtu', 'heat_2.vtu', 'heat_3.vtu', 'heat_4.vtu'], &
                             names=['piece-1', 'piece-2', 'piece-3', 'piece-4'], name='domain')
error = assembly%write_block(filenames=['probes.vtp'], names=['probes'], name='sensors')
error = assembly%write_block(action='close')
error = assembly%finalize()

Read back, the tree of blocks and datasets:

! read the assembly back: its blocks and datasets, depth first
error = assembly%initialize(filename='heat.vtm', action='read')
error = assembly%get_entries(level=level, kind=kind, name=name, file=file)
do e=1, size(level)
  print '(A,A,1X,A,1X,A)', repeat('  ', level(e) - 1), trim(kind(e)), trim(name(e)), trim(file(e))
enddo
error = assembly%finalize()

heat_7.f90

program heat
!< Tutorial, chapter 7: an assembly of the pieces of chapter 6 and of a set of probes, in a multi-block file.
use penf, only : I4P, R8P
use vtk_fortran, only : vtk_file, vtm_file
implicit none
integer(I4P), parameter :: probes=8
real(R8P),    parameter :: px(probes)=[0.2_R8P, 0.35_R8P, 0.5_R8P, 0.65_R8P, 0.8_R8P, 0.35_R8P, 0.7_R8P, 0.5_R8P]
real(R8P),    parameter :: py(probes)=[0.2_R8P, 0.2_R8P, 0.2_R8P, 0.2_R8P, 0.2_R8P, 0.4_R8P, 0.65_R8P, 0.8_R8P]
real(R8P),    parameter :: pz(probes)=0.5_R8P
real(R8P)               :: tp(probes)
type(vtk_file)          :: a_vtk_file
type(vtm_file)          :: assembly
integer(I4P),     allocatable :: level(:)
character(len=:), allocatable :: kind(:), name(:), file(:)
integer(I4P)            :: e, error

! the temperature at the probes, read from the pieces written by chapter 6
tp = probe(px, py, pz)

! the probes: a polydata of points, one vertex each
error = a_vtk_file%initialize(format='ascii', filename='probes.vtp', mesh_topology='PolyData')
error = a_vtk_file%xml_writer%write_piece(np=probes, nverts=probes, nlines=0, nstrips=0, npolys=0)
error = a_vtk_file%xml_writer%write_geo(np=probes, nc=probes, x=px, y=py, z=pz)
error = a_vtk_file%xml_writer%write_polydata_cells(verts_connectivity=[(e - 1, e=1, probes)], verts_offset=[(e, e=1, probes)])
error = a_vtk_file%xml_writer%write_dataarray(location='node', action='open')
error = a_vtk_file%xml_writer%write_dataarray(data_name='temperature', x=tp)
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 = assembly%initialize(filename='heat.vtm')
error = assembly%write_block(action='open', name='solver')
error = assembly%write_block(filenames=['heat_1.vtu', 'heat_2.vtu', 'heat_3.vtu', 'heat_4.vtu'], &
                             names=['piece-1', 'piece-2', 'piece-3', 'piece-4'], name='domain')
error = assembly%write_block(filenames=['probes.vtp'], names=['probes'], name='sensors')
error = assembly%write_block(action='close')
error = assembly%finalize()

! read the assembly back: its blocks and datasets, depth first
error = assembly%initialize(filename='heat.vtm', action='read')
error = assembly%get_entries(level=level, kind=kind, name=name, file=file)
do e=1, size(level)
  print '(A,A,1X,A,1X,A)', repeat('  ', level(e) - 1), trim(kind(e)), trim(name(e)), trim(file(e))
enddo
error = assembly%finalize()
contains
  impure elemental function probe(x, y, z) result(t)
  !< The temperature at the point nearest to (x,y,z), read from the piece that owns it.
  real(R8P), intent(in)  :: x, y, z
  real(R8P)              :: t
  real(R8P), allocatable :: qx(:), qy(:), qz(:), temperature(:)
  type(vtk_file)         :: piece
  character(len=10)      :: filename
  integer(I4P)           :: p, error

  ! the pieces split the cube along x in 4 equal parts: the owner of x
  p = min(4, int(x*4) + 1)
  write(filename, '(A,I1,A)') 'heat_', p, '.vtu'
  error = piece%initialize(filename=filename, action='read')
  error = piece%xml_reader%read_geo(x=qx, y=qy, z=qz)
  error = piece%xml_reader%read_dataarray(location='node', data_name='temperature', x=temperature)
  error = piece%finalize()
  t = temperature(minloc((qx - x)**2 + (qy - y)**2 + (qz - z)**2, dim=1))
  endfunction probe
endprogram heat

Running it

$ heat
block solver 
  block domain 
    dataset piece-1 heat_1.vtu
    dataset piece-2 heat_2.vtu
    dataset piece-3 heat_3.vtu
    dataset piece-4 heat_4.vtu
  block sensors 
    dataset probes probes.vtp

The assembly in ParaView, cut at the height of the probes:

What you learned

write_block with action='open'/'close' nests blocks; get_entries reads the tree back; a piece is read with action='read' and read_geo, read_dataarray. Reference: Nested blocks, Reading files.

Next: 8. Restart.