vtk_output Subroutine

  subroutine vtk_output (this, filename, iostat)
    use vtk_datasets,   only : unstruct_grid
    use vtk,            only : vtk_serial_write
    use vtk_cells,      only : voxel, vtkcell_list
    use vtk_attributes, only : attributes
    use array_functions_interface, only : OPERATOR(.catColumns.), OPERATOR(.columnVectors.)
    implicit none
    type vtk_data
      real(r8k), dimension(:), allocatable :: point_scalars, point_div_flux
    end type
    type(vtk_data), dimension(:), allocatable :: vtk_data_
    class(problem_discretization), intent(in) ::this
    character(LEN=*), intent(in), optional :: filename
    integer, intent(out) :: iostat
    type(voxel) :: voxel_cell  !! Voxel cell type
    type(vtkcell_list), dimension(:), allocatable :: cell_list !! Full list of all cells
    integer(i4k) :: ip, jp, kp !! block-local point id
    integer(i4k) :: ic, jc, kc !! block-local cell id
    real(r8k), dimension(:,:), allocatable :: points
    integer(i4k), dimension(:), allocatable :: block_cell_material, point_block_id
    type(attributes) :: grid_point_attributes, cell_attributes
    integer :: b, s, first_point_in_block, first_cell_in_block

    if (assertions) call assert(allocated(this%vertices), "problem_discretization%vtk_output: allocated(this%vertices())")

    allocate( cell_list(sum( this%vertices%num_cells())) )
    associate(my_first_block => lbound(this%vertices,1))
      allocate( points(this%vertices(my_first_block)%space_dimension(),0), block_cell_material(0), point_block_id(0) )
    end associate
    allocate( vtk_data_(this%num_scalars()) )
    do s = 1, this%num_scalars()
      allocate( vtk_data_(s)%point_scalars(0), vtk_data_(s)%point_div_flux(0) )
    end do

    first_point_in_block = 0
    first_cell_in_block = 1

    loop_over_grid_blocks: do b=lbound(this%vertices,1), ubound(this%vertices,1)

      associate( vertex_positions => this%vertices(b)%vectors() ) !! 4D array of nodal position vectors
        associate( npoints => shape(vertex_positions(:,:,:,1)) )  !! shape of 1st 3 dims specifies # points in ea. direction
          associate( ncells => npoints-1 )

            if (allocated(this%scalar_fields)) then
              loop_over_sclar_fields: &
              do s = 1, this%num_scalars()
                vtk_data_(s)%point_scalars = [vtk_data_(s)%point_scalars, this%scalar_fields(b,s)%get_scalar()]
              end do loop_over_sclar_fields
              if (allocated(this%scalar_flux_divergence)) then
                call assert(shape(this%scalar_fields)==shape(this%scalar_flux_divergence), "vtk_output: consistent scalar data")
                loop_over_flux_divergences: &
                do s = 1, this%num_scalars()
                  vtk_data_(s)%point_div_flux = [vtk_data_(s)%point_div_flux, this%scalar_flux_divergence(b,s)%get_scalar()]
                end do loop_over_flux_divergences
              end if
            end if

            points = points .catColumns. ( .columnVectors. vertex_positions )
            block_cell_material = [ block_cell_material, (this%vertices(b)%get_tag(), ic=1,PRODUCT(ncells)) ]
            point_block_id = [ point_block_id, [( b, ip=1, PRODUCT(npoints) )] ]

            do ic=1,ncells(1)
              do jc=1,ncells(2)
                do kc=1,ncells(3)
                  associate( block_local_point_id => & !! 8-element array of block-local IDs for voxel corners
                    [( [( [( kp*PRODUCT(npoints(1:2)) + jp*npoints(1) + ip, ip=ic,ic+1 )], jp=jc-1,jc )], kp=kc-1,kc )] &
                  )
                    call voxel_cell%setup ( first_point_in_block + block_local_point_id-1 )
                  end associate
                  associate( local_cell_id => (kc-1)*PRODUCT(ncells(1:2)) + (jc-1)*ncells(1) + ic )
                    associate(  cell_id => first_cell_in_block + local_cell_id - 1 )
                      allocate(cell_list(cell_id)%cell, source=voxel_cell)!GCC8 workaround for cell_list(i)%cell= voxel_cell
                    end associate
                  end associate
                end do
              end do
            end do

            first_cell_in_block  = first_cell_in_block  + PRODUCT( ncells  )
            first_point_in_block = first_point_in_block + PRODUCT( npoints )
          end associate
        end associate
      end associate
    end do loop_over_grid_blocks

    call assert( SIZE(point_block_id,1) == SIZE(points,2), "VTK block point data set & point set size match" )
    call assert( SIZE(block_cell_material,1) == SIZE(cell_list,1), "VTK cell data set & cell set size match" )

    call define_scalar(  cell_attributes, REAL( block_cell_material, r8k),  'material' )
    call define_scalar(  grid_point_attributes, REAL( point_block_id, r8k),  'block' )

    block
      type(unstruct_grid) vtk_grid
      type(attributes), dimension(:), allocatable :: scalar_attributes, flux_divergence_attributes, point_attributes
      character(len=*), parameter :: max_num_scalars='999'
      character(len=len(max_num_scalars)) scalar_number

      point_attributes = [grid_point_attributes]

      allocate( scalar_attributes(this%num_scalars()) )
      do s=1,size(scalar_attributes)
        write(scalar_number,'(i3)') s
        call define_scalar( scalar_attributes(s), vtk_data_(s)%point_scalars, 'scalar '//adjustl(trim(scalar_number)) )
        point_attributes = [point_attributes, scalar_attributes]
      end do

      allocate( flux_divergence_attributes(this%num_scalar_flux_divergences()) )
      do s=1,size(flux_divergence_attributes)
        write(scalar_number,'(i3)') s
        call define_scalar(flux_divergence_attributes(s), vtk_data_(s)%point_div_flux, 'divergence '//adjustl(trim(scalar_number)))
        point_attributes = [point_attributes, flux_divergence_attributes]
      end do

      call vtk_grid%init (points=points, cell_list=cell_list )
      call vtk_serial_write( &
        filename=filename, geometry=vtk_grid, multiple_io=.TRUE., celldatasets=[cell_attributes], pointdatasets=point_attributes)
    end block

    iostat = 0

#ifdef FORD
  end subroutine vtk_output