text: /******************************************************************************
   *
   * MantaFlow fluid solver framework
   * Copyright 2011-2016 Tobias Pfaff, Nils Thuerey  
   *
   * This program is free software, distributed under the terms of the
   * Apache License, Version 2.0 
   * http://www.apache.org/licenses/LICENSE-2.0
   *
   * Loading and writing grids and meshes to disk
   *
   ******************************************************************************/
  
  #include <iostream>
  #include <fstream>
  #include <cstdlib>
  #include <cstring>
  
  #if NO_ZLIB!=1
  extern "C" { 
  #include <zlib.h>
  }
  #endif
  
  #if OPENVDB==1
  #include "openvdb/openvdb.h"
  #endif
  
  #include "cnpy.h"
  #include "mantaio.h"
  #include "grid.h"
  #include "vector4d.h"
  #include "grid4d.h"
  
  using namespace std;
  
  namespace Manta {
  
  static const int STR_LEN_GRID  = 252;
  
  //! uni file header, v4
  typedef struct {
  	int dimX, dimY, dimZ; // grid size
  	int gridType, elementType, bytesPerElement; // data type info
  	char info[STR_LEN_GRID]; // mantaflow build information
  	int dimT;                // optionally store forth dimension for 4d grids
  	unsigned long long timestamp; // creation time
  } UniHeader;
  
  // note: header v4 only uses 4 bytes of the info string to store the fourth dimension, not needed for pdata
  
  
  //*****************************************************************************
  // conversion functions for double precision
  // (note - uni files always store single prec. values)
  //*****************************************************************************
  
  #if NO_ZLIB!=1
  template <class GRIDT> 
  void gridConvertWrite(gzFile& gzf, GRIDT& grid, void* ptr, UniHeader& head) {
  	errMsg("gridConvertWrite: unknown type, not yet supported");
  }
  
  template <>
  void gridConvertWrite(gzFile& gzf, Grid<int>& grid, void* ptr, UniHeader& head) {
  	gzwrite(gzf, &head,    sizeof(UniHeader));
  	gzwrite(gzf, &grid[0], sizeof(int)*head.dimX*head.dimY*head.dimZ);
  } 
  template <>
  void gridConvertWrite(gzFile& gzf, Grid<double>& grid, void* ptr, UniHeader& head) {
  	head.bytesPerElement = sizeof(float);
  	gzwrite(gzf, &head, sizeof(UniHeader));
  	float* ptrf = (float*)ptr;
  	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i,++ptrf) {
  		*ptrf = (float)grid[i];
  	} 
  	gzwrite(gzf, ptr, sizeof(float)* head.dimX*head.dimY*head.dimZ);
  } 
  template <>
  void gridConvertWrite(gzFile& gzf, Grid<Vector3D<double> >& grid, void* ptr, UniHeader& head) {
  	head.bytesPerElement = sizeof(Vector3D<float>);
  	gzwrite(gzf, &head, sizeof(UniHeader));
  	float* ptrf = (float*)ptr;
  	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i) {
  		for(int c=0; c<3; ++c) { *ptrf = (float)grid[i][c]; ptrf++; }
  	} 
  	gzwrite(gzf, ptr, sizeof(Vector3D<float>) *head.dimX*head.dimY*head.dimZ);
  }
  
  template <>
  void gridConvertWrite(gzFile& gzf, Grid4d<int>& grid, void* ptr, UniHeader& head) {
  	gzwrite(gzf, &head,    sizeof(UniHeader));
  	gzwrite(gzf, &grid[0], sizeof(int)*head.dimX*head.dimY*head.dimZ*head.dimT);
  }
  template <>
  void gridConvertWrite(gzFile& gzf, Grid4d<double>& grid, void* ptr, UniHeader& head) {
  	head.bytesPerElement = sizeof(float);
  	gzwrite(gzf, &head, sizeof(UniHeader));
  	float* ptrf = (float*)ptr;
  	IndexInt s = grid.getStrideT()*grid.getSizeT();
  	for(IndexInt i=0; i<s; ++i,++ptrf) {
  		*ptrf = (float)grid[i];
  	} 
  	gzwrite(gzf, ptr, sizeof(float)* s );
  } 
  template <>
  void gridConvertWrite(gzFile& gzf, Grid4d<Vector3D<double> >& grid, void* ptr, UniHeader& head) {
  	head.bytesPerElement = sizeof(Vector3D<float>);
  	gzwrite(gzf, &head, sizeof(UniHeader));
  	float* ptrf = (float*)ptr;
  	IndexInt s = grid.getStrideT()*grid.getSizeT();
  	for(IndexInt i=0; i<s; ++i) {
  		for(int c=0; c<3; ++c) { *ptrf = (float)grid[i][c]; ptrf++; }
  	} 
  	gzwrite(gzf, ptr, sizeof(Vector3D<float>) *s);
  }
  template <>
  void gridConvertWrite(gzFile& gzf, Grid4d<Vector4D<double> >& grid, void* ptr, UniHeader& head) {
  	head.bytesPerElement = sizeof(Vector4D<float>);
  	gzwrite(gzf, &head, sizeof(UniHeader));
  	float* ptrf = (float*)ptr;
  	IndexInt s = grid.getStrideT()*grid.getSizeT();
  	for(IndexInt i=0; i<s; ++i) {
  		for(int c=0; c<4; ++c) { *ptrf = (float)grid[i][c]; ptrf++; }
  	} 
  	gzwrite(gzf, ptr, sizeof(Vector4D<float>) *s);
  }
  
  template <class T>
  void gridReadConvert(gzFile& gzf, Grid<T>& grid, void* ptr, int bytesPerElement) {
  	errMsg("gridReadConvert: unknown type, not yet supported");
  }
  
  template <>
  void gridReadConvert<int>(gzFile& gzf, Grid<int>& grid, void* ptr, int bytesPerElement) {
  	gzread(gzf, ptr, sizeof(int)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ());
  	assertMsg (bytesPerElement == sizeof(int), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(int) );
  	// easy, nothing to do for ints
  	memcpy(&(grid[0]), ptr, sizeof(int) * grid.getSizeX()*grid.getSizeY()*grid.getSizeZ() );
  }
  
  template <>
  void gridReadConvert<double>(gzFile& gzf, Grid<double>& grid, void* ptr, int bytesPerElement) {
  	gzread(gzf, ptr, sizeof(float)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ());
  	assertMsg (bytesPerElement == sizeof(float), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(float) );
  	float* ptrf = (float*)ptr;
  	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i,++ptrf) {
  		grid[i] = (double)(*ptrf);
  	}
  }
  
  template <>
  void gridReadConvert<Vec3>(gzFile& gzf, Grid<Vec3>& grid, void* ptr, int bytesPerElement) {
  	gzread(gzf, ptr, sizeof(Vector3D<float>)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ());
  	assertMsg (bytesPerElement == sizeof(Vector3D<float>), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(Vector3D<float>) );
  	float* ptrf = (float*)ptr;
  	for(int i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i) {
  		Vec3 v;
  		for(int c=0; c<3; ++c) { v[c] = double(*ptrf); ptrf++; }
  		grid[i] = v;
  	}
  }
  
  template <class T>
  void gridReadConvert4d(gzFile& gzf, Grid4d<T>& grid, void* ptr, int bytesPerElement, int t) {
  	errMsg("gridReadConvert4d: unknown type, not yet supported");
  }
  
  template <>
  void gridReadConvert4d<int>(gzFile& gzf, Grid4d<int>& grid, void* ptr, int bytesPerElement, int t) { 
  	gzread(gzf, ptr, sizeof(int)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ() );
  	assertMsg (bytesPerElement == sizeof(int), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(int) );
  	// nothing to do for ints
  	memcpy(&(grid[grid.getSizeX()*grid.getSizeY()*grid.getSizeZ() *t]), ptr, sizeof(int) * grid.getSizeX()*grid.getSizeY()*grid.getSizeZ() );
  }
  
  template <>
  void gridReadConvert4d<double>(gzFile& gzf, Grid4d<double>& grid, void* ptr, int bytesPerElement, int t) {
  	assertMsg (bytesPerElement == sizeof(float), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(float) );
  
  	float* ptrf = (float*)ptr;
  	gzread(gzf, ptr, sizeof(float)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ());
  	for(IndexInt i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i , ++ptrf ) {
  		grid[ grid.getSizeX()*grid.getSizeY()*grid.getSizeZ()*t + i] = (double)(*ptrf);
  	} 
  }
  
  template <>
  void gridReadConvert4d<Vec3>(gzFile& gzf, Grid4d<Vec3>& grid, void* ptr, int bytesPerElement, int t) {
  	assertMsg (bytesPerElement == sizeof(Vector3D<float>), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(float) );
  
  	gzread(gzf, ptr, sizeof(Vector3D<float>)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ());
  	float* ptrf = (float*)ptr;
  	for(IndexInt i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i) {
  		Vec3 v;
  		for(int c=0; c<3; ++c) { v[c] = double(*ptrf); ptrf++; }
  		grid[grid.getSizeX()*grid.getSizeY()*grid.getSizeZ()*t + i] = v;
  	} 
  }
  
  template <>
  void gridReadConvert4d<Vec4>(gzFile& gzf, Grid4d<Vec4>& grid, void* ptr, int bytesPerElement, int t) {
  	assertMsg (bytesPerElement == sizeof(Vector4D<float>), "grid element size doesn't match "<< bytesPerElement <<" vs "<< sizeof(float) );
  
  	gzread(gzf, ptr, sizeof(Vector4D<float>)*grid.getSizeX()*grid.getSizeY()*grid.getSizeZ());
  	float* ptrf = (float*)ptr;
  	for(IndexInt i=0; i<grid.getSizeX()*grid.getSizeY()*grid.getSizeZ(); ++i) {
  		Vec4 v;
  		for(int c=0; c<4; ++c) { v[c] = double(*ptrf); ptrf++; }
  		grid[grid.getSizeX()*grid.getSizeY()*grid.getSizeZ()*t + i] = v;
  	} 
  }
  
  
  
  // make sure compatible grid types dont lead to errors...
  static int unifyGridType(int type) {
  	// real <> levelset
  	if(type & GridBase::TypeReal)     type |= GridBase::TypeLevelset;
  	if(type & GridBase::TypeLevelset) type |= GridBase::TypeReal;
  	// vec3 <> mac
  	if(type & GridBase::TypeVec3)     type |= GridBase::TypeMAC;
  	if(type & GridBase::TypeMAC)      type |= GridBase::TypeVec3;
  	return type;
  }
  
  #endif // NO_ZLIB!=1
  
  //*****************************************************************************
  // grid data
  //*****************************************************************************
  
  template<class T>
  void writeGridTxt(const string& name, Grid<T>* grid) {
  	debMsg( "writing grid " << grid->getName() << " to text file " << name ,1);
  
  	ofstream ofs(name.c_str());
  	if (!ofs.good())
  		errMsg("writeGridTxt: can't open file " << name);
  	FOR_IJK(*grid) {
  		ofs << Vec3i(i,j,k) <<" = "<< (*grid)(i,j,k) <<"\n";
  	}
  	ofs.close();
  }
  
  template<class T>
  void writeGridRaw(const string& name, Grid<T>* grid) {
  	debMsg( "writing grid " << grid->getName() << " to raw file " << name ,1);
  	
  #	if NO_ZLIB!=1
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "wb1"); // do some compression
  	if (!gzf) errMsg("writeGridRaw: can't open file " << name);
  	gzwrite(gzf, &((*grid)[0]), sizeof(T)*grid->getSizeX()*grid->getSizeY()*grid->getSizeZ());
  	gzclose(gzf);
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  }
  
  template<class T>
  void readGridRaw(const string& name, Grid<T>* grid) {
  	debMsg( "reading grid " << grid->getName() << " from raw file " << name ,1);
  	
  #	if NO_ZLIB!=1
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "rb");
  	if (!gzf) errMsg("readGridRaw: can't open file " << name);
  	
  	IndexInt bytes = sizeof(T)*grid->getSizeX()*grid->getSizeY()*grid->getSizeZ();
  	IndexInt readBytes = gzread(gzf, &((*grid)[0]), bytes);
  	assertMsg(bytes==readBytes, "can't read raw file, stream length does not match, "<<bytes<<" vs "<<readBytes);
  	gzclose(gzf);
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  }
  
  //! legacy headers for reading old files
  typedef struct {
  	int dimX, dimY, dimZ;
  	int frames, elements, elementType, bytesPerElement, bytesPerFrame;
  } UniLegacyHeader;
  
  typedef struct {
  	int dimX, dimY, dimZ;
  	int gridType, elementType, bytesPerElement;
  } UniLegacyHeader2;
  
  typedef struct {
  	int dimX, dimY, dimZ;
  	int gridType, elementType, bytesPerElement;
  	char info[256];
  	unsigned long long timestamp;
  } UniLegacyHeader3;
  
  //! for auto-init & check of results of test runs , optionally returns info string of header
  void getUniFileSize(const string& name, int& x, int& y, int& z, int* t, std::string* info) {
  	x = y = z = 0;
  #	if NO_ZLIB!=1
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "rb");
  	if (gzf) { 
  		char ID[5]={0,0,0,0,0};
  		gzread(gzf, ID, 4); 
  
  		// v3
  		if ( (!strcmp(ID, "MNT2")) || (!strcmp(ID, "M4T2")) ) {
  			UniLegacyHeader3 head;
  			assertMsg (gzread(gzf, &head, sizeof(UniLegacyHeader3)) == sizeof(UniLegacyHeader3), "can't read file, no header present"); 
  			x = head.dimX;
  			y = head.dimY;
  			z = head.dimZ;
  
  			// optionally , read fourth dim
  			if ((!strcmp(ID, "M4T2")) && t) {
  				int dimT = 0;
  				gzread(gzf, &dimT, sizeof(int) );
  				(*t) = dimT;
  			}
  		}
  
  		// v4
  		if ( (!strcmp(ID, "MNT3")) || (!strcmp(ID, "M4T3")) ) {
  			UniHeader head;
  			assertMsg (gzread(gzf, &head, sizeof(UniHeader)) == sizeof(UniHeader), "can't read file, no header present"); 
  			x = head.dimX;
  			y = head.dimY;
  			z = head.dimZ;
  			if(t) (*t) = head.dimT;
  		}
  
  		gzclose(gzf);
  	}
  #	endif
  	if(info) {
  		std::ostringstream out;
  		out << x <<","<< y <<","<< z ;
  		if(t && (*t)>0 ) out << ","<< (*t);
  		*info = out.str();
  	}
  }
  PYTHON() Vec3 getUniFileSize(const string& name) {
  	int x,y,z;
  	getUniFileSize(name, x,y,z);
  	return Vec3( Real(x), Real(y), Real(z) );
  }
  
  //! for test run debugging
  PYTHON() void printUniFileInfoString(const string& name) {
  	std::string info("<file not found>");
  	int x=-1,y=-1,z=-1,t=-1;
  	// use getUniFileSize to parse the different headers
  	getUniFileSize(name, x,y,z,&t, &info);
  	debMsg("File '"<<name<<"' info: "<< info ,1);
  }
  
  
  // actual read/write functions
  
  template <class T>
  void writeGridUni(const string& name, Grid<T>* grid) {
  	debMsg( "Writing grid " << grid->getName() << " to uni file " << name ,1);
  	
  #	if NO_ZLIB!=1
  	char ID[5] = "MNT3";
  	UniHeader head;
  	head.dimX = grid->getSizeX();
  	head.dimY = grid->getSizeY();
  	head.dimZ = grid->getSizeZ();
  	head.dimT = 0;
  	head.gridType = grid->getType();
  	head.bytesPerElement = sizeof(T);
  	snprintf( head.info, STR_LEN_GRID, "%s", buildInfoString().c_str() );	
  	MuTime stamp;
  	head.timestamp = stamp.time;
  	
  	if (grid->getType() & GridBase::TypeInt)
  		head.elementType = 0;
  	else if (grid->getType() & GridBase::TypeReal)
  		head.elementType = 1;
  	else if (grid->getType() & GridBase::TypeVec3)
  		head.elementType = 2;
  	else 
  		errMsg("writeGridUni: unknown element type");
  	
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "wb1"); // do some compression
  	if (!gzf) errMsg("writeGridUni: can't open file " << name);
  	
  	gzwrite(gzf, ID, 4);
  #	if FLOATINGPOINT_PRECISION!=1
  	// always write float values, even if compiled with double precision...
  	Grid<T> temp(grid->getParent());
  	// "misuse" temp grid as storage for floating point values (we have double, so it will always fit)
  	gridConvertWrite( gzf, *grid, &(temp[0]), head);
  #	else
  	void* ptr = &((*grid)[0]);
  	gzwrite(gzf, &head, sizeof(UniHeader));
  	gzwrite(gzf, ptr, sizeof(T)*head.dimX*head.dimY*head.dimZ);
  #	endif
  	gzclose(gzf);
  
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  };
  
  
  template <class T>
  void readGridUni(const string& name, Grid<T>* grid) {
  	debMsg( "Reading grid " << grid->getName() << " from uni file " << name ,1);
  
  #	if NO_ZLIB!=1
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "rb");
  	if (!gzf) errMsg("readGridUni: can't open file " << name);
  
  	char ID[5]={0,0,0,0,0};
  	gzread(gzf, ID, 4);
  	
  	if (!strcmp(ID, "DDF2")) {
  		// legacy file format
  		UniLegacyHeader head;
  		assertMsg (gzread(gzf, &head, sizeof(UniLegacyHeader)) == sizeof(UniLegacyHeader), "can't read file, no header present");
  		assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match");
  		assertMsg (head.bytesPerElement * head.elements == sizeof(T), "grid type doesn't match");
  		// skip flags
  		int numEl = head.dimX*head.dimY*head.dimZ;
  		gzseek(gzf, numEl, SEEK_CUR);
  		// actual grid read
  		gzread(gzf, &((*grid)[0]), sizeof(T)*numEl);
  	} 
  	else if (!strcmp(ID, "MNT1")) {
  		// legacy file format 2
  		UniLegacyHeader2 head;
  		assertMsg (gzread(gzf, &head, sizeof(UniLegacyHeader2)) == sizeof(UniLegacyHeader2), "can't read file, no header present");
  		assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< grid->getSize() );
  		assertMsg (head.gridType == grid->getType(), "grid type doesn't match "<< head.gridType<<" vs "<< grid->getType() );
  		assertMsg (head.bytesPerElement == sizeof(T), "grid element size doesn't match "<< head.bytesPerElement <<" vs "<< sizeof(T) );
  		gzread(gzf, &((*grid)[0]), sizeof(T)*head.dimX*head.dimY*head.dimZ);
  	}
  	else if (!strcmp(ID, "MNT2")) {
  		// a bit ugly, almost identical to MNT3
  		UniLegacyHeader3 head;
  		assertMsg (gzread(gzf, &head, sizeof(UniLegacyHeader3)) == sizeof(UniLegacyHeader3), "can't read file, no header present");
  		assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< grid->getSize() );
  		assertMsg (unifyGridType(head.gridType)==unifyGridType(grid->getType()) , "grid type doesn't match "<< head.gridType<<" vs "<< grid->getType() );
  #		if FLOATINGPOINT_PRECISION!=1
  		Grid<T> temp(grid->getParent());
  		void*  ptr  = &(temp[0]);
  		gridReadConvert<T>(gzf, *grid, ptr, head.bytesPerElement);
  #		else
  		assertMsg (head.bytesPerElement == sizeof(T), "grid element size doesn't match "<< head.bytesPerElement <<" vs "<< sizeof(T) );
  		gzread(gzf, &((*grid)[0]), sizeof(T)*head.dimX*head.dimY*head.dimZ);
  #		endif
  	} 
  	else if (!strcmp(ID, "MNT3")) {
  		// current file format
  		UniHeader head;
  		assertMsg (gzread(gzf, &head, sizeof(UniHeader)) == sizeof(UniHeader), "can't read file, no header present");
  		assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< grid->getSize() );
  		assertMsg (unifyGridType(head.gridType)==unifyGridType(grid->getType()) , "grid type doesn't match "<< head.gridType<<" vs "<< grid->getType() );
  #		if FLOATINGPOINT_PRECISION!=1
  		// convert float to double
  		Grid<T> temp(grid->getParent());
  		void*  ptr  = &(temp[0]);
  		gridReadConvert<T>(gzf, *grid, ptr, head.bytesPerElement);
  #		else
  		assertMsg (head.bytesPerElement == sizeof(T), "grid element size doesn't match "<< head.bytesPerElement <<" vs "<< sizeof(T) );
  		gzread(gzf, &((*grid)[0]), sizeof(T)*head.dimX*head.dimY*head.dimZ);
  #		endif
  	} else {
  		errMsg( "readGridUni: Unknown header '"<<ID<<"' " );
  	}
  	gzclose(gzf);
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  };
  
  template <class T>
  void writeGridVol(const string& name, Grid<T>* grid) {
  	debMsg( "writing grid " << grid->getName() << " to vol file " << name ,1);
  	errMsg("writeGridVol: Type not yet supported!");
  }
  
  struct volHeader { 
  	char 	ID[3];
  	char 	version;
  	int		encoding;
  	int		dimX, dimY, dimZ;
  	int		channels;
  	Vec3	bboxMin, bboxMax;
  };
  
  template <>
  void writeGridVol<Real>(const string& name, Grid<Real>* grid) {
  	debMsg( "writing real grid " << grid->getName() << " to vol file " << name ,1);
  	
  	volHeader header;
  	header.ID[0] = 'V';
  	header.ID[1] = 'O';
  	header.ID[2] = 'L';
  	header.version  = 3;
  	header.encoding = 1; // float32 precision
  	header.dimX = grid->getSizeX();
  	header.dimY = grid->getSizeY();
  	header.dimZ = grid->getSizeZ();
  	header.channels = 1; // only 1 channel
  	header.bboxMin = Vec3(-0.5);
  	header.bboxMax = Vec3( 0.5);
  
  	FILE* fp = fopen( name.c_str(), "wb" );
  	if (fp == NULL) {
  		errMsg("writeGridVol: Cannot open '" << name << "'");
  		return;
  	}
  
  	fwrite( &header, sizeof(volHeader), 1, fp );
  
  #	if FLOATINGPOINT_PRECISION==1
  	// for float, write one big chunk
  	fwrite( &(*grid)[0], sizeof(float), grid->getSizeX()*grid->getSizeY()*grid->getSizeZ(), fp );
  #	else
  	// explicitly convert each entry to float - we might have double precision in mantaflow
  	FOR_IDX(*grid) {
  		float value = (*grid)[idx];
  		fwrite( &value, sizeof(float), 1, fp );
  	}
  #	endif
  
  	fclose(fp);
  };
  
  
  template <class T>
  void readGridVol(const string& name, Grid<T>* grid) {
  	debMsg( "writing grid " << grid->getName() << " to vol file " << name ,1);
  	errMsg("readGridVol: Type not yet supported!");
  }
  
  template <>
  void readGridVol<Real>(const string& name, Grid<Real>* grid) {
  	debMsg( "reading real grid " << grid->getName() << " from vol file " << name ,1);
  	
  	volHeader header; 
  	FILE* fp = fopen( name.c_str(), "rb" );
  	if (fp == NULL) {
  		errMsg("readGridVol: Cannot open '" << name << "'");
  		return;
  	}
  
  	// note, only very basic file format checks here!
  	assertMsg( fread( &header, 1, sizeof(volHeader) ,  fp) == sizeof(volHeader), "can't read file, no header present");
  	if( header.dimX != grid->getSizeX() || header.dimY != grid->getSizeY() || header.dimZ != grid->getSizeZ()) errMsg( "grid dim doesn't match, "<< Vec3(header.dimX,header.dimY,header.dimZ)<<" vs "<< grid->getSize() );
  #	if FLOATINGPOINT_PRECISION!=1
  	errMsg("readGridVol: Double precision not yet supported");
  #	else
  	const unsigned int s = sizeof(float)*header.dimX*header.dimY*header.dimZ;
  	assertMsg( fread( &((*grid)[0]), 1, s, fp) == s, "can't read file, no / not enough data");
  #	endif
  
  	fclose(fp);
  };
  
  // 4d grids IO
  
  template <class T>
  void writeGrid4dUni(const string& name, Grid4d<T>* grid) {
  	debMsg( "writing grid4d " << grid->getName() << " to uni file " << name ,1);
  	
  #	if NO_ZLIB!=1
  	char ID[5] = "M4T3";
  	UniHeader head;
  	head.dimX = grid->getSizeX();
  	head.dimY = grid->getSizeY();
  	head.dimZ = grid->getSizeZ();
  	head.dimT = grid->getSizeT();
  	head.gridType = grid->getType();
  	head.bytesPerElement = sizeof(T);
  	snprintf( head.info, STR_LEN_GRID, "%s", buildInfoString().c_str() );	
  	MuTime stamp; stamp.get();
  	head.timestamp = stamp.time;
  	
  	if (grid->getType() & Grid4dBase::TypeInt)
  		head.elementType = 0;
  	else if (grid->getType() & Grid4dBase::TypeReal)
  		head.elementType = 1;
  	else if (grid->getType() & Grid4dBase::TypeVec3)
  		head.elementType = 2;
  	else if (grid->getType() & Grid4dBase::TypeVec4)
  		head.elementType = 2;
  	else 
  		errMsg("writeGrid4dUni: unknown element type");
  	
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "wb1"); // do some compression
  	if (!gzf) errMsg("writeGrid4dUni: can't open file " << name);
  	
  	gzwrite(gzf, ID, 4);
  #	if FLOATINGPOINT_PRECISION!=1
  	Grid4d<T> temp(grid->getParent());
  	gridConvertWrite< Grid4d<T> >( gzf, *grid, &(temp[0]), head);
  #	else
  	gzwrite(gzf, &head, sizeof(UniHeader));
  
  	// can be too large - write in chunks
  	for(int t=0; t<head.dimT; ++t) { 
  		void* ptr = &((*grid)[           head.dimX*head.dimY*head.dimZ* t ]);
  		gzwrite(gzf, ptr,      sizeof(T)*head.dimX*head.dimY*head.dimZ* 1);
  	}
  #	endif
  	gzclose(gzf);
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  };
  
  //! note, reading 4d uni grids is slightly more complicated than 3d ones
  //! as it optionally supports sliced reading
  template <class T>
  void readGrid4dUni(const string& name, Grid4d<T>* grid, int readTslice, Grid4d<T>* slice, void** fileHandle ) 
  {
  	if(grid)  debMsg( "reading grid "  << grid->getName()  << " from uni file " << name ,1);
  	if(slice) debMsg( "reading slice " << slice->getName() << ",t="<<readTslice<<" from uni file " << name ,1);
  
  #	if NO_ZLIB!=1
  	gzFile gzf = NULL;
  	char ID[5]={0,0,0,0,0};
  
  	// optionally - reuse file handle, if valid one is passed in fileHandle pointer...
  	if( (!fileHandle) || (fileHandle && (*fileHandle == NULL)) ) {
  		gzf = (gzFile) safeGzopen(name.c_str(), "rb");
  		if (!gzf) errMsg("readGrid4dUni: can't open file "<<name);
  
  		gzread(gzf, ID, 4);
  		if( fileHandle) { *fileHandle = gzf; }
  	} else {
  		// optimized read - reduced sanity checks
  		gzf = (gzFile)(*fileHandle); 
  		void* ptr = &( (*slice)[ 0 ] );
  		gzread(gzf, ptr, sizeof(T)* slice->getStrideT()* 1);  // quick and dirty...
  		return;
  	}
  	
  	if( (!strcmp(ID, "M4T2")) || (!strcmp(ID, "M4T3")) ) {
  		int headerSize = -1;
  
  		// current file format
  		UniHeader head;
  		if(!strcmp(ID, "M4T3")) {
  			headerSize = sizeof(UniHeader);
  			assertMsg (gzread(gzf, &head, sizeof(UniHeader)) == sizeof(UniHeader), "can't read file, no 4d header present");
  			if(FLOATINGPOINT_PRECISION==1) assertMsg (head.bytesPerElement == sizeof(T), "4d grid element size doesn't match "<< head.bytesPerElement <<" vs "<< sizeof(T) );
  		}
  		// old header
  		if(!strcmp(ID, "M4T2")) {
  			UniLegacyHeader3 lhead;
  			headerSize = sizeof(UniLegacyHeader3) + sizeof(int);
  			assertMsg (gzread(gzf, &lhead, sizeof(UniLegacyHeader3)) == sizeof(UniLegacyHeader3), "can't read file, no 4dl header present");
  			if(FLOATINGPOINT_PRECISION==1) assertMsg (lhead.bytesPerElement == sizeof(T), "4d grid element size doesn't match "<< lhead.bytesPerElement <<" vs "<< sizeof(T) ); 
  
  			int fourthDim  = 0;
  			gzread(gzf, &fourthDim, sizeof(fourthDim));
  
  			head.dimX = lhead.dimX;
  			head.dimY = lhead.dimY;
  			head.dimZ = lhead.dimZ;
  			head.dimT = fourthDim;
  			head.gridType = lhead.gridType;
  		}
  
  		if(readTslice<0) {
  			assertMsg (head.dimX == grid->getSizeX() && head.dimY == grid->getSizeY() && head.dimZ == grid->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< grid->getSize() );
  			assertMsg ( unifyGridType(head.gridType)==unifyGridType(grid->getType()) , "grid type doesn't match "<< head.gridType<<" vs "<< grid->getType() );
  
  			// read full 4d grid
  			assertMsg (head.dimT == grid->getSizeT(), "grid dim4 doesn't match, "<< head.dimT <<" vs "<< grid->getSize() );
  
  			// can be too large - read in chunks
  #			if FLOATINGPOINT_PRECISION!=1
  			Grid4d<T> temp(grid->getParent());
  			void* ptr = &(temp[0]); 
  			for(int t=0; t<head.dimT; ++t) { 
  				gridReadConvert4d<T>(gzf, *grid, ptr, head.bytesPerElement, t);
  			}
  #			else
  			for(int t=0; t<head.dimT; ++t) { 
  				void* ptr = &((*grid)[            head.dimX*head.dimY*head.dimZ* t ] );
  				gzread(gzf, ptr,       sizeof(T)* head.dimX*head.dimY*head.dimZ*    1); 
  			}
  #			endif
  		} else {
  			// read chosen slice only
  			assertMsg (head.dimX == slice->getSizeX() && head.dimY == slice->getSizeY() && head.dimZ == slice->getSizeZ(), "grid dim doesn't match, "<< Vec3(head.dimX,head.dimY,head.dimZ)<<" vs "<< slice->getSize() );
  			assertMsg ( unifyGridType(head.gridType)==unifyGridType(slice->getType()) , "grid type doesn't match "<< head.gridType<<" vs "<< slice->getType() );
  
  #			if FLOATINGPOINT_PRECISION!=1
  			errMsg( "readGrid4dUni: NYI (2)" ); // slice read not yet supported for double
  #			else
  			assertMsg( slice, "No 3d slice grid data given" );
  			assertMsg (readTslice < head.dimT, "grid dim4 slice too large "<< readTslice <<" vs "<< head.dimT );
  			void* ptr = &( (*slice)[ 0 ] );
  			gzseek(gzf, sizeof(T)*head.dimX*head.dimY*head.dimZ* readTslice + headerSize + 4 , SEEK_SET );
  			gzread(gzf, ptr,         sizeof(T)*head.dimX*head.dimY*head.dimZ* 1);
  #			endif
  		}
  	} else {
  		debMsg( "Unknown header!" ,1);
  	}
  
  	if( !fileHandle) { 
  		gzclose(gzf);
  	}
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  };
  void readGrid4dUniCleanup(void** fileHandle) {
  	gzFile gzf = NULL; 
  	if( fileHandle) {
  		gzf = (gzFile)(*fileHandle);
  		gzclose(gzf);
  		*fileHandle = NULL;
  	}
  }
  
  template<class T>
  void writeGrid4dRaw(const string& name, Grid4d<T>* grid) {
  	debMsg( "writing grid4d " << grid->getName() << " to raw file " << name ,1);
  	
  #	if NO_ZLIB!=1
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "wb1"); // do some compression
  	if (!gzf) errMsg("writeGrid4dRaw: can't open file " << name);
  	gzwrite(gzf, &((*grid)[0]), sizeof(T)*grid->getSizeX()*grid->getSizeY()*grid->getSizeZ()*grid->getSizeT());
  	gzclose(gzf);
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  }
  
  template<class T>
  void readGrid4dRaw(const string& name, Grid4d<T>* grid) {
  	debMsg( "reading grid4d " << grid->getName() << " from raw file " << name ,1);
  	
  #	if NO_ZLIB!=1
  	gzFile gzf = (gzFile) safeGzopen(name.c_str(), "rb");
  	if (!gzf) errMsg("readGrid4dRaw: can't open file " << name);
  	
  	IndexInt bytes = sizeof(T)*grid->getSizeX()*grid->getSizeY()*grid->getSizeZ()*grid->getSizeT();
  	IndexInt readBytes = gzread(gzf, &((*grid)[0]), bytes);
  	assertMsg(bytes==readBytes, "can't read raw file, stream length does not match, "<<bytes<<" vs "<<readBytes);
  	gzclose(gzf);
  #	else
  	debMsg( "file format not supported without zlib" ,1);
  #	endif
  }
  
  
  //*****************************************************************************
  // optional openvdb export
  
  #if OPENVDB==1
  
  template <class T>
  void writeGridVDB(const string& name, Grid<T>* grid) {
  	debMsg("Writing grid " << grid->getName() << " to vdb file " << name << " not yet supported!", 1);
  }
  
  template <class T>
  void readGridVDB(const string& name, Grid<T>* grid) {
  	debMsg("Reading grid " << grid->getName() << " from vdb file " << name << " not yet supported!", 1);
  }
  
  template <>
  void writeGridVDB(const string& name, Grid<int>* grid) {
  	debMsg("Writing int grid " << grid->getName() << " to vdb file " << name, 1);
  
  	// Create an empty int32-point grid with background value 0.
  	openvdb::initialize();
  	openvdb::Int32Grid::Ptr gridVDB = openvdb::Int32Grid::create();
  	gridVDB->setTransform( openvdb::math::Transform::createLinearTransform( 1./grid->getSizeX() )); //voxel size
  
  	// Get an accessor for coordinate-based access to voxels.
  	openvdb::Int32Grid::Accessor accessor = gridVDB->getAccessor();
  
  	gridVDB->setGridClass(openvdb::GRID_UNKNOWN);
  
  	// Name the grid "density".
  	gridVDB->setName( grid->getName() );
  
  	openvdb::io::File file(name);
  
  	FOR_IJK(*grid) {
  		openvdb::Coord xyz(i, j, k);
  		accessor.setValue(xyz, (*grid)(i, j, k));
  	}
  
  	// Add the grid pointer to a container.
  	openvdb::GridPtrVec gridsVDB;
  	gridsVDB.push_back(gridVDB);
  
  	// Write out the contents of the container.
  	file.write(gridsVDB);
  	file.close();
  }
  
  template <>
  void readGridVDB(const string& name, Grid<int>* grid) {
  	debMsg("Reading int grid " << grid->getName() << " from vdb file " << name, 1);
  
  	openvdb::initialize();
  	openvdb::io::File file(name);
  	file.open();
  
  	openvdb::GridBase::Ptr baseGrid;
  	for (openvdb::io::File::NameIterator nameIter = file.beginName(); nameIter != file.endName(); ++nameIter)
  	{
  	#ifndef BLENDER
  		// Read in only the grid we are interested in.
  		if (nameIter.gridName() == grid->getName()) {
  			baseGrid = file.readGrid(nameIter.gridName());
  		} else {
  			debMsg("skipping grid " << nameIter.gridName(), 1);
  		}
  	#else
  		// For Blender, skip name check and pick first grid from loop
  		baseGrid = file.readGrid(nameIter.gridName());
  		break;
  	#endif
  	}
  	file.close();
  	openvdb::Int32Grid::Ptr gridVDB = openvdb::gridPtrCast<openvdb::Int32Grid>(baseGrid);
  
  	openvdb::Int32Grid::Accessor accessor = gridVDB->getAccessor();
  
  	FOR_IJK(*grid) {
  		openvdb::Coord xyz(i, j, k);
  		int v = accessor.getValue(xyz);
  		(*grid)(i, j, k) = v;
  	}
  }
  
  template <>
  void writeGridVDB(const string& name, Grid<Real>* grid) {
  	debMsg("Writing real grid " << grid->getName() << " to vdb file " << name, 1);
  
  	// Create an empty floating-point grid with background value 0.
  	openvdb::initialize();
  	openvdb::FloatGrid::Ptr gridVDB = openvdb::FloatGrid::create(); 
  	gridVDB->setTransform( openvdb::math::Transform::createLinearTransform( 1./grid->getSizeX() )); //voxel size
  
  	// Get an accessor for coordinate-based access to voxels.
  	openvdb::FloatGrid::Accessor accessor = gridVDB->getAccessor();
  
  	// Identify the grid as a level set.
  	gridVDB->setGridClass(openvdb::GRID_FOG_VOLUME);
  
  	// Name the grid "density".
  	gridVDB->setName( grid->getName() );
  
  	openvdb::io::File file(name);
  
  	FOR_IJK(*grid) { 
  		openvdb::Coord xyz(i, j, k);
  		accessor.setValue(xyz, (*grid)(i, j, k));	
  	}
  
  	// Add the grid pointer to a container.
  	openvdb::GridPtrVec gridsVDB;
  	gridsVDB.push_back(gridVDB);
  
  	// Write out the contents of the container.
  	file.write(gridsVDB);
  	file.close();
  };
  
  template <>
  void readGridVDB(const string& name, Grid<Real>* grid) {
  	debMsg("Reading real grid " << grid->getName() << " from vdb file " << name, 1);
  
  	openvdb::initialize();
  	openvdb::io::File file(name);
  	file.open();
  
  	openvdb::GridBase::Ptr baseGrid;
  	for (openvdb::io::File::NameIterator nameIter = file.beginName(); nameIter != file.endName(); ++nameIter)
  	{
  	#ifndef BLENDER
  		// Read in only the grid we are interested in.
  		if (nameIter.gridName() == grid->getName()) {
  			baseGrid = file.readGrid(nameIter.gridName());
  		} else {
  			debMsg("skipping grid " << nameIter.gridName(), 1);
  		}
  	#else
  		// For Blender, skip name check and pick first grid from loop
  		baseGrid = file.readGrid(nameIter.gridName());
  		break;
  	#endif
  	}
  	file.close();
  	openvdb::FloatGrid::Ptr gridVDB = openvdb::gridPtrCast<openvdb::FloatGrid>(baseGrid);
  
  	openvdb::FloatGrid::Accessor accessor = gridVDB->getAccessor();
  
  	FOR_IJK(*grid) {
  		openvdb::Coord xyz(i, j, k);
  		float v = accessor.getValue(xyz);
  		(*grid)(i, j, k) = v;
  	}
  };
  
  template <>
  void writeGridVDB(const string& name, Grid<Vec3>* grid) {
  	debMsg("Writing vec3 grid " << grid->getName() << " to vdb file " << name, 1);
  
  	openvdb::initialize(); 
  	openvdb::Vec3SGrid::Ptr gridVDB = openvdb::Vec3SGrid::create();
  	// note , warning - velocity content currently not scaled...
  	gridVDB->setTransform( openvdb::math::Transform::createLinearTransform( 1./grid->getSizeX() )); //voxel size 
  	openvdb::Vec3SGrid::Accessor accessor = gridVDB->getAccessor();
  
  	// MAC or regular vec grid?
  	if(grid->getType() & GridBase::TypeMAC) gridVDB->setGridClass(openvdb::GRID_STAGGERED);
  	else gridVDB->setGridClass(openvdb::GRID_UNKNOWN);
  
  	gridVDB->setName( grid->getName() );
  
  	openvdb::io::File file(name);
  	FOR_IJK(*grid) { 
  		openvdb::Coord xyz(i, j, k);
  		Vec3 v = (*grid)(i, j, k);
  		openvdb::Vec3f vo( (float)v[0] , (float)v[1] , (float)v[2] );
  		accessor.setValue(xyz, vo);
  	}
  
  	openvdb::GridPtrVec gridsVDB;
  	gridsVDB.push_back(gridVDB);
  
  	file.write(gridsVDB);
  	file.close();
  };
  
  template <>
  void readGridVDB(const string& name, Grid<Vec3>* grid) {
  	debMsg("Reading vec3 grid " << grid->getName() << " from vdb file " << name, 1);
  
  	openvdb::initialize();
  	openvdb::io::File file(name);
  	file.open();
  
  	openvdb::GridBase::Ptr baseGrid;
  	for (openvdb::io::File::NameIterator nameIter = file.beginName(); nameIter != file.endName(); ++nameIter)
  	{
  	#ifndef BLENDER
  		// Read in only the grid we are interested in.
  		if (nameIter.gridName() == grid->getName()) {
  			baseGrid = file.readGrid(nameIter.gridName());
  		} else {
  			debMsg("skipping grid " << nameIter.gridName(), 1);
  		}
  	#else
  		// For Blender, skip name check and pick first grid from loop
  		baseGrid = file.readGrid(nameIter.gridName());
  		break;
  	#endif
  	}
  	file.close();
  	openvdb::Vec3SGrid::Ptr gridVDB = openvdb::gridPtrCast<openvdb::Vec3SGrid>(baseGrid);
  
  	openvdb::Vec3SGrid::Accessor accessor = gridVDB->getAccessor();
  
  	FOR_IJK(*grid) {
  		openvdb::Coord xyz(i, j, k);
  		openvdb::Vec3f v = accessor.getValue(xyz);
  		(*grid)(i, j, k).x = (float)v[0];
  		(*grid)(i, j, k).y = (float)v[1];
  		(*grid)(i, j, k).z = (float)v[2];
  	}
  };
  
  #endif // OPENVDB==1
  
  
  //*****************************************************************************
  // npz file support (warning - read works, but write generates uncompressed npz; i.e. not recommended for large volumes)
  
  template <class T>
  void writeGridNumpy(const string& name, Grid<T>* grid) {
  #	if NO_ZLIB==1
  	debMsg( "file format not supported without zlib" ,1);
  	return;
  #	endif
  #	if FLOATINGPOINT_PRECISION!=1
  	errMsg("writeGridNumpy: Double precision not yet supported");
  #	endif
  
  	// find suffix to differentiate between npy <-> npz , TODO: check for actual "npy" string
  	std::string::size_type idx;
  	bool bUseNpz = false;
  	idx = name.rfind('.');
  	if(idx != std::string::npos) {
  		bUseNpz = name.substr(idx+1) == "npz";
  		debMsg( "Writing grid " << grid->getName() << " to npz file " << name ,1);
  	} else {
  		debMsg( "Writing grid " << grid->getName() << " to npy file " << name ,1);
  	}
  
  	// storage code
  	size_t uDim = 1;
  	if (grid->getType() & GridBase::TypeInt || grid->getType() & GridBase::TypeReal || grid->getType() & GridBase::TypeLevelset)
  		uDim = 1;
  	else if (grid->getType() & GridBase::TypeVec3 || grid->getType() & GridBase::TypeMAC)
  		uDim = 3;
  	else
  		errMsg("writeGridNumpy: unknown element type");
  
  	const std::vector<size_t> shape = {static_cast<size_t>(grid->getSizeZ()), static_cast<size_t>(grid->getSizeY()), static_cast<size_t>(grid->getSizeX()), uDim};
  
  	if(bUseNpz){
  		// note, the following generates a zip file without compression
  		if (grid->getType() & GridBase::TypeVec3 || grid->getType() & GridBase::TypeMAC) {
  			// cast to float* for export!
  			float* ptr = (float*)&((*grid)[0]);
  			cnpy::npz_save(name, "arr_0", ptr, shape, "w");
  		} else {
  			T* ptr = &((*grid)[0]);
  			cnpy::npz_save(name, "arr_0", ptr, shape, "w");
  		}
  	} else {
  		cnpy::npy_save(name, &grid[0], shape, "w");
  	}
  };
  
  template <class T>
  void readGridNumpy(const string& name, Grid<T>* grid) {
  #	if NO_ZLIB==1
  	debMsg( "file format not supported without zlib" ,1);
  	return;
  #	endif
  #	if FLOATINGPOINT_PRECISION!=1
  	errMsg("readGridNumpy: Double precision not yet supported");
  #	endif
  
  	// find suffix to differentiate between npy <-> npz
  	std::string::size_type idx;
  	bool bUseNpz = false;
  	idx = name.rfind('.');
  	if(idx != std::string::npos) {
  		bUseNpz = name.substr(idx+1) == "npz";
  		debMsg( "Reading grid " << grid->getName() << " as npz file " << name ,1);
  	} else {
  		debMsg( "Reading grid " << grid->getName() << " as npy file " << name ,1);
  	}
  
      cnpy::NpyArray gridArr;
      if(bUseNpz) {
  		cnpy::npz_t fNpz = cnpy::npz_load(name);
  		gridArr = fNpz["arr_0"];
      } else {
          gridArr = cnpy::npy_load(name);
      }
  
      // Check the file meta information
      assertMsg (gridArr.shape[2] == grid->getSizeX() && gridArr.shape[1] == grid->getSizeY() && gridArr.shape[0] == grid->getSizeZ(), "grid dim doesn't match, " << Vec3(gridArr.shape[2], gridArr.shape[1], gridArr.shape[0]) << " vs " << grid->getSize() );
      size_t uDim = 1;
      if (grid->getType() & GridBase::TypeInt || grid->getType() & GridBase::TypeReal || grid->getType() & GridBase::TypeLevelset)
          uDim = 1;
      else if (grid->getType() & GridBase::TypeVec3 || grid->getType() & GridBase::TypeMAC)
          uDim = 3;
      else
          errMsg("readGridNumpy: unknown element type");
      assertMsg (gridArr.shape[3] == uDim, "grid data dim doesn't match, " << gridArr.shape[3] << " vs " << uDim );
  
  	if (grid->getType() & GridBase::TypeVec3 || grid->getType() & GridBase::TypeMAC) {
  		// treated as float* for export , thus consider 3 elements
  	    assertMsg (3*gridArr.word_size == sizeof(T), "vec3 grid data size doesn't match, " << 3*gridArr.word_size << " vs " << sizeof(T));
  	} else {
  	    assertMsg (gridArr.word_size == sizeof(T), "grid data size doesn't match, " << gridArr.word_size << " vs " << sizeof(T));
  	}
  
      // copy back, TODO: beautify...
     	memcpy(&((*grid)[0]), gridArr.data<T>(), sizeof(T) * grid->getSizeX() * grid->getSizeY() * grid->getSizeZ() );
  };
  
  // adopted from getUniFileSize
  void getNpzFileSize(const string& name, int& x, int& y, int& z, int* t = NULL, std::string* info = NULL) {
  	x = y = z = 0;
  #	if NO_ZLIB!=1
  	debMsg( "file format not supported without zlib" ,1);
  	return;
  #	endif
  #	if FLOATINGPOINT_PRECISION!=1
  	errMsg("getNpzFileSize: Double precision not yet supported");
  #	endif
  	// find suffix to differentiate between npy <-> npz
      cnpy::NpyArray gridArr;
  	cnpy::npz_t fNpz = cnpy::npz_load(name);
  	gridArr = fNpz["arr_0"];
  
  	z = gridArr.shape[0];
  	y = gridArr.shape[1];
  	x = gridArr.shape[2];
  	if(t) (*t) = 0; // unused for now
  }
  PYTHON() Vec3 getNpzFileSize(const string& name) {
  	int x,y,z;
  	getNpzFileSize(name, x,y,z); 
  	return Vec3( Real(x), Real(y), Real(z) );
  }
  
  
  //*****************************************************************************
  // helper functions
  
  
  void quantizeReal(Real& v, const Real step) { 
  	int    q  = int(v / step + step*0.5);
  	double qd = q * (double)step;
  	v = (Real)qd;
  }
  KERNEL(idx) void knQuantize(Grid<Real>& grid, Real step) {
  	quantizeReal( grid(idx), step );
  } 
  PYTHON() void quantizeGrid(Grid<Real>& grid, Real step) { knQuantize(grid,step); }
  
  KERNEL(idx) void knQuantizeVec3(Grid<Vec3>& grid, Real step) {
  	for(int c=0; c<3; ++c) quantizeReal( grid(idx)[c], step );
  } 
  PYTHON() void quantizeGridVec3(Grid<Vec3>& grid, Real step) { knQuantizeVec3(grid,step); }
  
  
  
  // explicit instantiation
  template void writeGridRaw<int> (const string& name, Grid<int>*  grid);
  template void writeGridRaw<Real>(const string& name, Grid<Real>* grid);
  template void writeGridRaw<Vec3>(const string& name, Grid<Vec3>* grid);
  template void writeGridUni<int> (const string& name, Grid<int>*  grid);
  template void writeGridUni<Real>(const string& name, Grid<Real>* grid);
  template void writeGridUni<Vec3>(const string& name, Grid<Vec3>* grid);
  template void writeGridVol<int> (const string& name, Grid<int>*  grid);
  template void writeGridVol<Vec3>(const string& name, Grid<Vec3>* grid);
  template void writeGridTxt<int> (const string& name, Grid<int>*  grid);
  template void writeGridTxt<Real>(const string& name, Grid<Real>* grid);
  template void writeGridTxt<Vec3>(const string& name, Grid<Vec3>* grid);
  
  template void readGridRaw<int>  (const string& name, Grid<int>*  grid);
  template void readGridRaw<Real> (const string& name, Grid<Real>* grid);
  template void readGridRaw<Vec3> (const string& name, Grid<Vec3>* grid);
  template void readGridUni<int>  (const string& name, Grid<int>*  grid);
  template void readGridUni<Real> (const string& name, Grid<Real>* grid);
  template void readGridUni<Vec3> (const string& name, Grid<Vec3>* grid);
  template void readGridVol<int>  (const string& name, Grid<int>*  grid);
  template void readGridVol<Vec3> (const string& name, Grid<Vec3>* grid);
  
  template void readGrid4dUni<int>  (const string& name, Grid4d<int>*  grid, int readTslice, Grid4d<int>*  slice, void** fileHandle);
  template void readGrid4dUni<Real> (const string& name, Grid4d<Real>* grid, int readTslice, Grid4d<Real>* slice, void** fileHandle);
  template void readGrid4dUni<Vec3> (const string& name, Grid4d<Vec3>* grid, int readTslice, Grid4d<Vec3>* slice, void** fileHandle);
  template void readGrid4dUni<Vec4> (const string& name, Grid4d<Vec4>* grid, int readTslice, Grid4d<Vec4>* slice, void** fileHandle);
  template void writeGrid4dUni<int> (const string& name, Grid4d<int>*  grid);
  template void writeGrid4dUni<Real>(const string& name, Grid4d<Real>* grid);
  template void writeGrid4dUni<Vec3>(const string& name, Grid4d<Vec3>* grid);
  template void writeGrid4dUni<Vec4>(const string& name, Grid4d<Vec4>* grid);
  
  template void readGrid4dRaw<int>  (const string& name, Grid4d<int>*  grid);
  template void readGrid4dRaw<Real> (const string& name, Grid4d<Real>* grid);
  template void readGrid4dRaw<Vec3> (const string& name, Grid4d<Vec3>* grid);
  template void readGrid4dRaw<Vec4> (const string& name, Grid4d<Vec4>* grid);
  template void writeGrid4dRaw<int> (const string& name, Grid4d<int>*  grid);
  template void writeGrid4dRaw<Real>(const string& name, Grid4d<Real>* grid);
  template void writeGrid4dRaw<Vec3>(const string& name, Grid4d<Vec3>* grid);
  template void writeGrid4dRaw<Vec4>(const string& name, Grid4d<Vec4>* grid);
  
  template void writeGridNumpy<int> (const string& name, Grid<int>*  grid);
  template void writeGridNumpy<Real>(const string& name, Grid<Real>* grid);
  template void writeGridNumpy<Vec3>(const string& name, Grid<Vec3>* grid);
  template void readGridNumpy<int>  (const string& name, Grid<int>*  grid);
  template void readGridNumpy<Real> (const string& name, Grid<Real>* grid);
  template void readGridNumpy<Vec3> (const string& name, Grid<Vec3>* grid);
  
  #if OPENVDB==1
  template void writeGridVDB<int>(const string& name, Grid<int>*  grid);
  template void writeGridVDB<Vec3>(const string& name, Grid<Vec3>* grid);
  template void writeGridVDB<Real>(const string& name, Grid<Real>* grid);
  
  template void readGridVDB<int>(const string& name, Grid<int>*  grid);
  template void readGridVDB<Vec3>(const string& name, Grid<Vec3>* grid);
  template void readGridVDB<Real>(const string& name, Grid<Real>* grid);
  #endif // OPENVDB==1
  
  } //namespace
  
  namespace Manta {
  
  }
  