#!/bin/bash

# This utility converts AHFinderDirect h.t*.ah*.gp files into h*_it*.vtk
# files (the iteration/horizon index rearrangement makes VisIt correctly
# identify the databases). 
# Thanks to Roland Haas for contributing the gawk loop.

if [ $# -ne 1 ]
then
  echo "Usage: gp2vtk <h.t*.ah*.gp file>"
else
  dir=`dirname $1`
  datafile=`basename $1`

  base=${datafile%*.ah*.gp}
  sid=${datafile%*.gp}
  sid=${sid#*ah}
  iid=${base#h.t*}
  vtkfile=h${sid}_it${iid}.vtk

  echo "Horizon ${sid}, iteration $iid"

  egrep -v '^$' $dir/$datafile | egrep -v '^#' | awk '{ print $4 "\t" $5 "\t" $6 }' > tmp
  N=`wc -l tmp | awk '{ print $1 }'`

  echo "# vtk DataFile Version 2.0" > $vtkfile
  echo "# Horizon data" >> $vtkfile
  echo "ASCII" >> $vtkfile
  echo "DATASET POLYDATA" >> $vtkfile
  echo "POINTS $N float" >> $vtkfile
  cat tmp >> $vtkfile
  \rm -rf tmp

  npoints=`cat $dir/$datafile | grep N_rho | head -1 | awk '{ print $4 }'`
  nblocks=6
  nlines=`echo "$N $nblocks / $npoints / p" | dc`
  maxi=`echo "$nlines 2 - p" | dc`
  maxj=`echo "$npoints 2 - p" | dc`
  npolygons=`echo "$nblocks $maxi 1 + * $maxj 1 + * p" | dc`
  nsize=`echo "$npolygons 5 * p" | dc`
  echo "" >> $vtkfile
  echo "POLYGONS $npolygons $nsize" >> $vtkfile

  gawk -vnblocks=$nblocks -vmaxi=$maxi -vmaxj=$maxj \
       -vnpoints=$npoints -vvtkfile=$vtkfile -vnlines=$nlines '
  BEGIN {
    OFS="\t"
    k=0
    while ( k < nblocks )
    {
  	i=0
  	while ( i <= maxi )
        {
  		j=0
  		printf("%s", "BLOCK "k" of "nblocks": "i"/"maxi" \r")
  		while ( j <= maxj )
  		{
  			v1=(k * nlines + i) * npoints + j
  			v2=v1 + 1
  			v4=v1 +  npoints
  			v3=v4 + 1
  			print 4,v1,v2,v3,v4 >>vtkfile
  			j=j+1
                  }
  		i=i+1
        }
  	print ""
  	k=k+1
    }
  }'

  echo "" >> $vtkfile
  echo "POINT_DATA $N" >> $vtkfile
fi