#!/usr/bin/perl
use strict;
use MIME::Base64;
use Compress::Raw::Zlib;

# Copyright (C) 2008 Stefan Tibus
#
# inp2vtu_gz.pl is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# inp2vtu_gz.pl is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with magpar; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

# 2007-11-20 STi: inp2vtu_gz.pl
# Supports node data and writes compressed binary data Vtk-XML files.
# Nodes and cells have to be ordered and continuous.
#
# 2008-02-01 STi: inp2vtu_gz.pl
# Added support for cell data.

# Note: Additional information on writing binary vtk files taken from
# http://news.tiker.net/node/430
# "What they don't tell you about VTK XML binary formats"

# Check Number of Arguments
# =========================
if($#ARGV!=0) {
	die("Usage: inp2vtu_gz.pl <problem>\n");
}

# Get Arguments
# =============
my $problem=$ARGV[0];

# Get List of .inp Files
# ======================
my @inpfiles=`ls -1 \"$problem.\"????\".inp\"`;
# Strip leading/trailing spaces/tabs/newlines
foreach my $inpfile (@inpfiles) {
	$inpfile=~s/^\s*(.*)\s*\n?/$1/;
}
# Skip last file if it has the number 9999
if(@inpfiles[-1]=~/^$problem\.9999\.inp/) {
	pop(@inpfiles);
}

#DBG#foreach my $inpfile (@inpfiles) {
#DBG#	print(stderr $inpfile."\n");
#DBG#}

# =================
# Create .vtu Files
# =================
my @vtufiles;

foreach my $inpfile (@inpfiles) {
	my $vtufile=$inpfile;
	# Change and store filename
	$vtufile=~s/(.*)\.(\d\d\d\d)\.inp/$1_$2.vtu/;
	push(@vtufiles,$vtufile);
	
	# skip conversion if VTU file is newer!
	if ( (stat($vtufile))[9] > (stat($inpfile))[9] ) {
	  print "Skipping $inpfile because $vtufile is newer!\n\n";
	  next;
	}

	# Read From .inp File
	# ===================
	open(INPFILE,"<",$inpfile) || die("Cannot open file \"$inpfile\" for read: $!\n");
	
	# Read Header
	# -----------
	# nnodes, ncells, nnodedata, ncelldata, nmodeldata
	my $line=<INPFILE>;
	$line=~s/^\s*(.*)\s*$/$1/;
	my($nnodes,$ncells,$nnodedata,$ncelldata,$nmodeldata)=split(/\s+/,$line,5);
	if(($nmodeldata!=0)) {
		die("Don't know how to handle model data in file \"$inpfile\"!\n");
	}

	# Read Nodes
	# ----------
	print(STDERR $nnodes," nodes\n");
	my @nodes;
	for(my $i=0;$i<$nnodes;$i++) {
		$line=<INPFILE>;
		$line=~s/^\s*(.*)\s*$/$1/;
		my($n,$x,$y,$z)=split(/\s+/,$line,4);
		if($n!=($i+1)) {
			die("Don't know how to handle unordered nodes in file \"$inpfile\"!\n");
		}
		push(@nodes,[$x,$y,$z]);
	}

	# Read Cells
	# ----------
	print(STDERR $ncells," cells\n");
	my @cells;
	for(my $i=0;$i<$ncells;$i++) {
		$line=<INPFILE>;
		$line=~s/^\s*(.*)\s*$/$1/;
		my($n,$mat,$type,$remainder)=split(/\s+/,$line,4);
		if($type!='tet') {
			die("Don't know how to handle cells other than tetrahedra in file \"$inpfile\"!\n");
		}
		if($n!=($i+1)) {
			die("Don't know how to handle unordered cells in file \"$inpfile\"!\n");
		}
		my($n0,$n1,$n2,$n3)=split(/\s+/,$remainder,4);
		push(@cells,[4,$n0,$n1,$n2,$n3,$mat]);
	}

	# Read Node Data
	# --------------
	my @nodelabels;
	my @nodedimensions;
	my @nodefirstelems;
	my @nodelabeldatatypes;
	my $nnodelabels;
	my @nodevectors;
    if($nnodedata>0) {
	print(STDERR $nnodedata," nodedata\n");
	$line=<INPFILE>;
	$line=~s/^\s*(.*)\s*$/$1/;
	my($nnodevalues,$remainder)=split(/\s+/,$line,2);
	my(@snodevalues)=split(/\s+/,$remainder,$nnodevalues);
	my $totalsize;
	for(my $i=0;$i<$nnodevalues;$i++) {
		$totalsize+=$snodevalues[$i];
	}
	if($totalsize!=$nnodedata) {
		die("Node data size mismatch in file \"$inpfile\"!\n");
	}
	# Read data labels
	my $label;
	my $dimension;
	my $firstelem=0;
	my $labeldatatype;
	for(my $i=0;$i<$nnodevalues;$i++) {
		$line=<INPFILE>;
		$line=~s/^\s*(.*)\s*$/$1/;
		$line=~s/^(.*),.*$/$1/;
		if($line=~/_?[XxYyZz]$/) {
			# expect x,y,z
			$line=~s/^(.*?)_?[XxYyZz]$/$1/;
			if($line eq $label) {
				# increase dimensionality of former label
				$nodedimensions[-1]++;
				# increase index of first element of next label
				$firstelem++;
				next;
			} else {
				$label=$line;
				$dimension=1;
			}
		} else {
			$label=$line;
			$dimension=1;
		}
		push(@nodelabels,$label);
		push(@nodedimensions,$dimension);
		push(@nodefirstelems,$firstelem);
		if($label eq "id") {
		    $labeldatatype="UInt32";
		} elsif($label eq "proc") {
		    $labeldatatype="UInt32";
		} elsif($label eq "vert_pid") {
		    $labeldatatype="UInt32";
		} elsif($label eq "prop") {
		    $labeldatatype="UInt32";
		} else {
		    $labeldatatype="Float32";
		}
		push(@nodelabeldatatypes,$labeldatatype);
		# increase index of first element of next label
		$firstelem++;
	}
	$nnodelabels=scalar(@nodelabels);
	# Read data values
	for(my $i=0;$i<$nnodes;$i++) {
		$line=<INPFILE>;
		$line=~s/^\s*(.*)\s*$/$1/;
		my($n,@vector)=split(/\s+/,$line,($nnodevalues+1));
		if($n!=($i+1)) {
			die("Don't know how to handle unordered node data in file \"$inpfile\"!\n");
		}
		push(@nodevectors,[@vector]);
	}
    }

	# Read Cell Data
	# --------------
	my @celllabels;
	my @celldimensions;
	my @cellfirstelems;
	my @celllabeldatatypes;
	my $ncelllabels;
	my @cellvectors;
    if($ncelldata>0) {
	print(STDERR $ncelldata," celldata\n");
	$line=<INPFILE>;
	$line=~s/^\s*(.*)\s*$/$1/;
	my($ncellvalues,$remainder)=split(/\s+/,$line,2);
	my(@scellvalues)=split(/\s+/,$remainder,$ncellvalues);
	my $totalsize;
	for(my $i=0;$i<$ncellvalues;$i++) {
		$totalsize+=$scellvalues[$i];
	}
	if($totalsize!=$ncelldata) {
		die("Cell data size mismatch in file \"$inpfile\"!\n");
	}
	# Read data labels
	my $label;
	my $dimension;
	my $firstelem=0;
	my $labeldatatype;
	for(my $i=0;$i<$ncellvalues;$i++) {
		$line=<INPFILE>;
		$line=~s/^\s*(.*)\s*$/$1/;
		$line=~s/^(.*),.*$/$1/;
		if($line=~/_?[XxYyZz]$/) {
			# expect x,y,z
			$line=~s/^(.*?)_?[XxYyZz]$/$1/;
			if($line eq $label) {
				# increase dimensionality of former label
				$celldimensions[-1]++;
				# increase index of first element of next label
				$firstelem++;
				next;
			} else {
				$label=$line;
				$dimension=1;
			}
		} else {
			$label=$line;
			$dimension=1;
		}
		push(@celllabels,$label);
		push(@celldimensions,$dimension);
		push(@cellfirstelems,$firstelem);
		if($label eq "id") {
		    $labeldatatype="UInt32";
		} elsif($label eq "proc") {
		    $labeldatatype="UInt32";
		} elsif($label eq "vert_pid") {
		    $labeldatatype="UInt32";
		} elsif($label eq "prop") {
		    $labeldatatype="UInt32";
		} else {
		    $labeldatatype="Float32";
		}
		push(@celllabeldatatypes,$labeldatatype);
		# increase index of first element of next label
		$firstelem++;
	}
	$ncelllabels=scalar(@celllabels);
	# Read data values
	for(my $i=0;$i<$ncells;$i++) {
		$line=<INPFILE>;
		$line=~s/^\s*(.*)\s*$/$1/;
		my($n,@vector)=split(/\s+/,$line,($ncellvalues+1));
		if($n!=($i+1)) {
			die("Don't know how to handle unordered cell data in file \"$inpfile\"!\n");
		}
		push(@cellvectors,[@vector]);
	}
    }

	# Close .inp File
	# ---------------
	close(INPFILE);

	# Helper Functions
	# ================

my %packformats=( "UInt8" => "C", "UInt16" => "S", "UInt32" => "L", "Float32" => "f" );
my %binsizes=( "UInt8" => 1, "UInt16" => 2, "UInt32" => 4, "Float32" => 4 );

	# Encode Binary Data
	# ------------------
	# Prepend appropriate header and base64-encode
sub encodebin {
	my $bindata=shift(@_);

	my $binsize=length($bindata);
	my $base64header=encode_base64(pack("L",$binsize));
	my $base64data=encode_base64($bindata);
	my $base64record=$base64header.$base64data;
	$base64record=~s/\n//g;
	$base64record.="\n";

	return($base64record);
}

	# Deflate Using Zlib
	# ------------------
sub deflatezlib {
	my $indata=shift(@_);
	my $status;
	my $outpart;
	my $outdata;

	my $deflator=new Compress::Raw::Zlib::Deflate or die("Cannot create a deflation stream!\n");
	$status=$deflator->deflate($indata,$outpart);
	if($status!=Z_OK) {
		die("Deflation failed!\n");
	}
	$outdata.=$outpart;
	$status=$deflator->flush($outpart);
	if($status!=Z_OK) {
		die("Deflation failed!\n");
	}
	$outdata.=$outpart;

	return($outdata);
}

	# Inflate Using Zlib
	# ------------------
sub inflatezlib {
	my $indata=shift(@_);
	my $status;
	my $outdata;

	my $inflator=new Compress::Raw::Zlib::Inflate or die("Cannot create an inflation stream!\n");
	$status=$inflator->inflate($indata,$outdata);
	if($status!=Z_STREAM_END) {
		die("Inflation failed!\n");
	}

	return($outdata);
}

sub encodebinzlib {
	my $bindata=shift(@_);

	my $binsize=length($bindata);
#DBG#print($binsize);
	my $zippeddata=deflatezlib($bindata);
	my $zippedsize=length($zippeddata);
#DBG#print(" -> ".$zippedsize);
	my $base64header=encode_base64(pack("L4",1,$binsize,$binsize,$zippedsize));
	my $base64data=encode_base64($zippeddata);
#DBG#print(" -> ".length($base64data));
#DBG#my $unbase64data=decode_base64($base64data);
#DBG#my $unzippeddata=inflatezlib($unbase64data);
#DBG#$base64header=encode_base64(pack("L",$binsize));
#DBG#$base64data=encode_base64($unzippeddata);
#DBG#print(" -> ".length($unzippeddata)."\n");
	my $base64record=$base64header.$base64data;
	$base64record=~s/\n//g;
	$base64record.="\n";

	return($base64record);
}

	# Verify Data Type Sizes
	# ======================
	if(length(pack("c",1))!=1) {
		die("char is not 8 bits long!\n");
	}
	if(length(pack("C",1))!=1) {
		die("unsigned char is not 8 bits long!\n");
	}
	if(length(pack("s",1))!=2) {
		die("short is not 16 bits long!\n");
	}
	if(length(pack("S",1))!=2) {
		die("unsigned short is not 16 bits long!\n");
	}
	if(length(pack("l",1))!=4) {
		die("long is not 32 bits long!\n");
	}
	if(length(pack("L",1))!=4) {
		die("unsigned long is not 32 bits long!\n");
	}
	if(length(pack("f",1.0))!=4) {
		die("float is not 32 bits long!\n");
	}

	# Write to .vtu File
	# ==================
	open(VTUFILE,">",$vtufile) || die("Cannot open or create file \"$vtufile\" for write: $!\n");
	binmode(VTUFILE,':raw');
	# Write Header
	# ------------
	print(VTUFILE "<?xml version=\"1.0\"?>\n");
	print(VTUFILE "<VTKFile type=\"UnstructuredGrid\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n");
	# Note: Intel 80x86 is little endian, Sun SPARC and IBM PowerPC are big endian platforms.
	# See: http://www.netrino.com/Publications/Glossary/Endianness.php
	# See: http://www.intel.com/design/intarch/papers/endian.pdf
	print(VTUFILE "  <UnstructuredGrid>\n");
	print(VTUFILE "    <Piece NumberOfPoints=\"$nnodes\" NumberOfCells=\"$ncells\">\n");

	# Write Point (Node) Data
	# -----------------------
	# User data
	# Set $nnodelabels=1 if you want first data set only (data for mumag / M for magpar)
	#$nnodelabels=1;
	if($nnodelabels>0) {
		if($nodedimensions[0]>1) {
			print(VTUFILE "      <PointData Vectors=\"$nodelabels[0]\">\n");
		} else {
			print(VTUFILE "      <PointData>\n");
		}
			for(my $i=0;$i<$nnodelabels;$i++) {
				my $label=$nodelabels[$i];
				my $dimension=$nodedimensions[$i];
				my $firstelem=$nodefirstelems[$i];
				my $labeldatatype=$nodelabeldatatypes[$i];
				my $labelpackformat=$packformats{$labeldatatype};
				my $bindata;
				if($dimension==1) {
					print(VTUFILE "        <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"binary\">\n");
					for(my $j=0;$j<$nnodes;$j++) {
						$bindata.=pack( $labelpackformat,      $nodevectors[$j][$firstelem]);	
					}
				} elsif($dimension==2) {
					print(VTUFILE "        <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"binary\">\n");
					for(my $j=0;$j<$nnodes;$j++) {
						$bindata.=pack( $labelpackformat."2",           $nodevectors[$j][$firstelem],$nodevectors[$j][$firstelem+1]);	
					}
				} elsif($dimension==3) {
					print(VTUFILE "        <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"binary\">\n");
					for(my $j=0;$j<$nnodes;$j++) {
						$bindata.=pack( $labelpackformat."3",                 $nodevectors[$j][$firstelem],$nodevectors[$j][$firstelem+1],$nodevectors[$j][$firstelem+2]);	
					}
				} else {
					die("Don't know how to handle vectors of dimension > 3 in file \"$inpfile\"!\n");
				}
				my $binsize=length($bindata);
				if($binsize!=($dimension*$nnodes*$binsizes{$labeldatatype})) {
					die("Size of binary data for \"$label\" is wrong!\n");
				}
				print(VTUFILE encodebinzlib($bindata));
				print(VTUFILE "        </DataArray>\n");
	}
		print(VTUFILE "      </PointData>\n");
	}

	# Write Cell Data
	# ---------------
	# Material
	print(VTUFILE "      <CellData Scalars=\"Material\">\n");
	print(VTUFILE "        <DataArray Name=\"Material\" NumberOfComponents=\"1\" type=\"UInt32\" format=\"binary\">\n");
	my $bindata;
	for (my $i=0;$i<$ncells;$i++) {
		$bindata.=pack( $packformats{"UInt32"},  $cells[$i][5]);	
	}
	my $binsize=length($bindata);
	if($binsize!=$ncells*$binsizes{"UInt32"}) {
		die("Size of binary data for \"Material\" is wrong!\n");
	}
	print(VTUFILE encodebinzlib($bindata));
	print(VTUFILE "        </DataArray>\n");
	# User data
	# Set $ncelllabels=0 if you want to suppress cell data
	#$ncelllabels=0;
	if($ncelllabels>0) {
		for(my $i=0;$i<$ncelllabels;$i++) {
			my $label=$celllabels[$i];
			my $dimension=$celldimensions[$i];
			my $firstelem=$cellfirstelems[$i];
			my $labeldatatype=$celllabeldatatypes[$i];
			my $labelpackformat=$packformats{$labeldatatype};
			my $bindata;
			if($dimension==1) {
				print(VTUFILE "        <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"binary\">\n");
				for(my $j=0;$j<$ncells;$j++) {
					$bindata.=pack( $labelpackformat,      $cellvectors[$j][$firstelem]);	
				}
			} elsif($dimension==2) {
				print(VTUFILE "        <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"binary\">\n");
				for(my $j=0;$j<$ncells;$j++) {
					$bindata.=pack( $labelpackformat."2",           $cellvectors[$j][$firstelem],$cellvectors[$j][$firstelem+1]);	
				}
			} elsif($dimension==3) {
				print(VTUFILE "        <DataArray Name=\"$label\" NumberOfComponents=\"$dimension\" type=\"$labeldatatype\" format=\"binary\">\n");
				for(my $j=0;$j<$ncells;$j++) {
					$bindata.=pack( $labelpackformat."3",                 $cellvectors[$j][$firstelem],$cellvectors[$j][$firstelem+1],$cellvectors[$j][$firstelem+2]);	
				}
			} else {
				die("Don't know how to handle vectors of dimension > 3 in file \"$inpfile\"!\n");
			}
			my $binsize=length($bindata);
			if($binsize!=($dimension*$ncells*$binsizes{$labeldatatype})) {
				die("Size of binary data for \"$label\" is wrong!\n");
			}
			print(VTUFILE encodebinzlib($bindata));
			print(VTUFILE "        </DataArray>\n");
		}
	}
	print(VTUFILE "      </CellData>\n");

	# Write Nodes
	# -----------
	# Point coordinates
	print(VTUFILE "      <Points>\n");
	print(VTUFILE "        <DataArray NumberOfComponents=\"3\" type=\"Float32\" format=\"binary\">\n");
	my $bindata;
	for(my $i=0;$i<$nnodes;$i++) {
		$bindata.=pack( "f3",                 $nodes[$i][0],$nodes[$i][1],$nodes[$i][2]);	
	}
	my $binsize=length($bindata);
	if($binsize!=$nnodes*3*4) {
		die("Size of binary data for points is wrong!\n");
	}
	print(VTUFILE encodebinzlib($bindata));
	print(VTUFILE "        </DataArray>\n");
	print(VTUFILE "      </Points>\n");

	# Write Cells
	# -----------
	# Cell connectivity (tetrahedrons only)
	print(VTUFILE "      <Cells>\n");
	print(VTUFILE "        <DataArray Name=\"connectivity\" type=\"UInt32\" format=\"binary\">\n");
	my $bindata;
	for(my $i=0;$i<$ncells;$i++) {
		$bindata.=pack( "L4",               $cells[$i][1]-1,$cells[$i][2]-1,$cells[$i][3]-1,$cells[$i][4]-1);
	}
	my $binsize=length($bindata);
	if($binsize!=$ncells*4*4) {
		die("Size of binary data for \"connectivity\" is wrong!\n");
	}
	print(VTUFILE encodebinzlib($bindata));
	print(VTUFILE "        </DataArray>\n");
	# Cell connectivity offsets (for tetrahedrons)
	print(VTUFILE "        <DataArray Name=\"offsets\" type=\"UInt32\" format=\"binary\">\n");
	my $bindata;
	for(my $i=0;$i<$ncells;$i++) {
		$bindata.=pack( "L",  ($i+1)*4);
	}
	my $binsize=length($bindata);
	if($binsize!=$ncells*4) {
		die("Size of binary data for \"offsets\" is wrong!\n");
	}
	print(VTUFILE encodebinzlib($bindata));
	print(VTUFILE "        </DataArray>\n");
	# Cell types (tetrahedrons only)
	# Tetrahedron = 10
	print(VTUFILE "        <DataArray Name=\"types\" type=\"UInt8\" format=\"binary\">\n");
	my $bindata;
	for(my $i=0;$i<$ncells;$i++) {
		$bindata.=pack( "C",  10);
	}
	my $binsize=length($bindata);
	if($binsize!=$ncells*1) {
		die("Size of binary data for \"types\" is wrong!\n");
	}
	print(VTUFILE encodebinzlib($bindata));
	print(VTUFILE "        </DataArray>\n");
	print(VTUFILE "      </Cells>\n");

	# Close .vtu File
	# ---------------
	print(VTUFILE "    </Piece>\n");
	print(VTUFILE "  </UnstructuredGrid>\n");
	print(VTUFILE "</VTKFile>\n");
	close(VTUFILE);
	print($vtufile."\n");
}

# ================
# Create .pvd File
# ================

if(scalar(@vtufiles)>1) {
    my $pvdfile=$problem.".pvd";

    open(PVDFILE,">",$pvdfile) || die("Cannot open or create file \"$pvdfile\" for write: $!\n");

    # Write Header
    # ============
    print(PVDFILE "<?xml version=\"1.0\"?>\n");
    print(PVDFILE "<VTKFile type=\"Collection\">\n");
    print(PVDFILE "  <Collection>\n");

    # Get List of .vtu Files
    # ======================
    my @vtufiles=`ls -1 \"${problem}_\"????\".vtu\"`;
    # Strip leading/trailing spaces/tabs/newlines
    foreach my $vtufile (@vtufiles) {
	    $vtufile=~s/^\s*(.*)\s*\n?/$1/;
    }

    # Write List of .vtu Files
    # ========================
    my $i=0;
    foreach my $vtufile (@vtufiles) {
	    printf(PVDFILE "    <DataSet timestep=\"%04d\" file=\"%s\"/>\n",$i,$vtufile);
	    $i++;
    }

    # Close .pvd File
    # ===============
    print(PVDFILE "  </Collection>\n");
    print(PVDFILE "</VTKFile>\n");
    close(PVDFILE);
    print($pvdfile."\n");
}

# vim:set ft=perl sts=4 sw=4:
