mesh_vtk_file.cpp

#ifdef _MSC_VER //kill some warnings
#define _CRT_SECURE_NO_WARNINGS
#endif


#include "inmost.h"
#include <cfloat>

#if defined(USE_MESH)

//vtk states

#define R_VERSION     0
#define R_USERDATA    1
#define R_DATATYPE    2
#define R_WAITDATA    3

#define R_SGRID       4
#define R_RGRID       5
#define R_UGRID       6
#define R_POLYDATA    7
#define R_FIELD       8
#define R_SPOINTS     9

#define R_ATTRIBUTES  10
#define R_ATTRDATA    11

#define R_QUIT		  100

//static int __isnan__(double x) { return x != x; }
//static int isinf(double x) { return !isnan(x) && isnan(x - x); }
//static int __isinf__(double x) { return fabs(x) > DBL_MAX; }
//static int __isbad(double x) { return __isnan__(x) || __isinf__(x); }



template<typename T>
void ReadCoords(FILE * f,INMOST_DATA_REAL_TYPE c[3])
{
	T temp[3];
	if( fread(temp,sizeof(T),3,f) == 3 )
	{
		c[0] = temp[0];
		c[1] = temp[1];
		c[2] = temp[2];
	}
	else 
	{
		std::cout << __FILE__ << ":" << __LINE__ << " Cannot read 3 coordinates from file " << std::endl;
		throw INMOST::BadFile;
	}
}

namespace INMOST
{


	void skip_empty(char * readline,FILE * f)
	{
		bool empty = true;
		do
		{
			for(size_t k = 0; k < strlen(readline); ++k)
				if( !isspace(readline[k]) ) empty = false;
			if( empty )
			{
				if( fgets(readline,2048,f) == NULL ) 
					break;
			}
		} while(empty);
		
	}

	void skip_metadata(int comps, char * readline,FILE * f)
	{
		int dcnt,filled, scnt;
		//char ntype[256];
		if( !strncmp(readline,"METADATA",8) )
		{
			for(int k = 0; k < 2; ++k)
			{
				if( fgets(readline,2048,f) == NULL ) throw BadFile;
				if( !strncmp(readline,"INFORMATION",11) )
				{
					filled = sscanf(readline,"%*s %d",&dcnt); //INFORMATION %d
					if( filled != 1 ) throw BadFile;
					for(int k = 0; k < dcnt; ++k) //2 lines per data
					{
						if( fgets(readline,2048,f) == NULL )  //NAME
							throw BadFile;
						filled = fscanf(f,"%*s"); //DATA
						filled = fscanf(f,"%d",&scnt); //vector size
						if( filled != 1 ) throw BadFile;
						for(int q = 0; q < scnt; ++q) filled = fscanf(f,"%*s\n");
					}
				}
				else if( !strncmp(readline,"COMPONENT_NAMES",15) )
				{
					for(int k = 0; k < comps; ++k) //2 lines per data
					{
						if( fgets(readline,2048,f) == NULL ) 
							throw BadFile;
					}
				}
			}
		}
	}

	std::string Space2Underscore(const std::string & inp)
	{
		std::string ret = inp;
		for (size_t k = 0; k < ret.size(); ++k)
		{
			if (ret[k] == ' ') ret[k] = '_';
		}
		return ret;
	}
  
	int VtkElementType(ElementType t)
	{
		switch(t)
		{
			case Element::Line: return 3;
			case Element::Tri: return 5;
			case Element::Quad: return 9;
			case Element::MultiLine: return 4;
			case Element::Polygon: return 7;
			case Element::Tet: return 10;
			case Element::Hex: return 12;
			case Element::Prism: return 13;
			case Element::Pyramid: return 14;
			case Element::Polyhedron: return 42;
			case Element::MultiPolygon: return 42;
		}
		assert(false);
		return -1;
	}
	
	INMOST_DATA_ENUM_TYPE VtkElementNodes(ElementType t)
	{
		switch(t)
		{
			case Element::Line: return 2;
			case Element::Tri: return 3;
			case Element::Quad: return 4;
			case Element::MultiLine: return ENUMUNDEF;
			case Element::Polygon: return ENUMUNDEF;
			case Element::Tet: return 4;
			case Element::Hex: return 8;
			case Element::Prism: return 6;
			case Element::Pyramid: return 5;
			case Element::Polyhedron: return ENUMUNDEF;
			case Element::MultiPolygon: return ENUMUNDEF;
		}
		assert(false);
		return ENUMUNDEF;
	}
	

  void Mesh::SaveVTK(std::string File)
  {
		integer dim = GetDimensions();
		if( dim > 3 )
		{
			std::cout << "VTK file supports 3 dimensions max" << std::endl;
			return;
		}
		std::ofstream f(File.c_str());
		//~ FILE * f = fopen(File.c_str(),"w");
		//~ if( !f ) throw BadFileName;
		if( f.fail() ) throw BadFileName;
		f << "# vtk DataFile Version 3.0\n";
		if (this->GetFileOption("VTK_OUTPUT_FACES") == "1") f << "VTK_OUTPUT_FACES file is written by INMOST\n";
		else f << "file is written by INMOST\n";
		f << "ASCII\n";
		f << "DATASET UNSTRUCTURED_GRID\n";

		bool output_faces = false;
		for (INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
		{
			if (file_options[k].first == "VTK_OUTPUT_FACES")
			{
				output_faces = true;
			}
		}
		
		std::set< std::string > nosave, saveonly;
		
		nosave = TagOptions("nosave");
		saveonly = TagOptions("saveonly");
		
		//ReorderEmpty(CELL | NODE);
		Tag set_id = CreateTag("TEMPORARY_ELEMENT_ID",DATA_INTEGER,CELL |FACE| NODE,NONE,1);
		Storage::integer cur_num = 0;
		for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it) it->IntegerDF(set_id) = cur_num++;
		cur_num = 0;
		for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); ++it) it->IntegerDF(set_id) = cur_num++;
		cur_num = 0;

		if (output_faces) for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); ++it) it->IntegerDF(set_id) = cur_num++;

		f << "POINTS " << NumberOfNodes() << " double\n";
		for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); it++)
		{
			Storage::real_array coords = it->RealArray(CoordsTag());
			for(integer i = 0; i < dim; i++) 
			{
				double temp = coords[i];
				f << temp << " ";
			}
			for(integer i = dim; i < 3; i++)
				f << "0 ";
			f << "\n";
		}
		{
			std::vector<int> values;
			for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
			{	
				switch(it->GetGeometricType())
				{
					case Element::Tri:
					case Element::Quad:
					//case Element::MultiLine:
					//case Element::Polygon:
					case Element::Tet:
					case Element::Hex:
					case Element::Pyramid:
					case Element::Prism:
					{
						ElementArray<Node> nodes = it->getNodes();
						if( nodes.size() != VtkElementNodes(it->GetGeometricType()) ) goto safe_output;
						values.push_back(static_cast<integer>(nodes.size()));
						for(ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
							values.push_back(jt->IntegerDF(set_id));
						break;
					}
					case Element::MultiLine:
					case Element::Polygon:
					{
						ElementArray<Node> nodes = it->getNodes();
						values.push_back(static_cast<integer>(nodes.size()));
						for(ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
							values.push_back(jt->IntegerDF(set_id));
						break;
					}
						/*
					case Element::Prism:
					{
						ElementArray<Node> nodes = it->getNodes();
						if( nodes.size() != 6 ) goto safe_output;
						values.push_back(static_cast<integer>(nodes.size()));
						values.push_back(nodes[0].IntegerDF(set_id));
						values.push_back(nodes[2].IntegerDF(set_id));
						values.push_back(nodes[1].IntegerDF(set_id));
						values.push_back(nodes[3].IntegerDF(set_id));
						values.push_back(nodes[5].IntegerDF(set_id));
						values.push_back(nodes[4].IntegerDF(set_id));
						break;
					}
						 */
					/*
					case Element::Hex:
					{
						ElementArray<Node> nodes = it->getNodes();
						if( nodes.size() != 8 ) goto safe_output;
						values.push_back(static_cast<integer>(nodes.size()));
						values.push_back(nodes[0].IntegerDF(set_id));
						values.push_back(nodes[3].IntegerDF(set_id));
						values.push_back(nodes[2].IntegerDF(set_id));
						values.push_back(nodes[1].IntegerDF(set_id));
						values.push_back(nodes[4].IntegerDF(set_id));
						values.push_back(nodes[7].IntegerDF(set_id));
						values.push_back(nodes[6].IntegerDF(set_id));
						values.push_back(nodes[5].IntegerDF(set_id));
						break;
					}
					case Element::Pyramid:
					{
						ElementArray<Node> nodes = it->getNodes();
						if( nodes.size() != 5 ) goto safe_output;
						values.push_back(static_cast<integer>(nodes.size()));
						values.push_back(nodes[0].IntegerDF(set_id));
						values.push_back(nodes[3].IntegerDF(set_id));
						values.push_back(nodes[2].IntegerDF(set_id));
						values.push_back(nodes[1].IntegerDF(set_id));
						values.push_back(nodes[4].IntegerDF(set_id));
						break;
					}
					 */
					case Element::Polyhedron:
					case Element::MultiPolygon:
					{
safe_output:
						
						ElementArray<Face> faces = it->getFaces();
						integer totalNum = 1 + static_cast<integer>(faces.size());
                        for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); jt++)
							totalNum += jt->nbAdjElements(NODE);
						values.push_back(totalNum);
						values.push_back(static_cast<integer>(faces.size()));
						for(ElementArray<Face>::iterator jt = faces.begin(); jt != faces.end(); jt++)
						{
							ElementArray<Node> nodes = jt->getNodes();
							values.push_back(static_cast<integer>(nodes.size()));
							if( jt->FaceOrientedOutside(it->self()) )
								for(ElementArray<Node>::iterator kt = nodes.begin(); kt != nodes.end(); kt++)
									values.push_back(kt->IntegerDF(set_id));
							else
								for(ElementArray<Node>::reverse_iterator kt = nodes.rbegin(); kt != nodes.rend(); kt++)
									values.push_back(kt->IntegerDF(set_id));
						}
						break;
					}
					default: std::cout << __FILE__ << ":" << __LINE__ << " This should not happen " << Element::GeometricTypeName(it->GetGeometricType()) << std::endl;
				}
			}
			if (output_faces)
			{
				for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); it++)
				{
					switch (it->GetGeometricType())
					{
					case Element::Line:
					{
						ElementArray<Node> nodes = it->getNodes();
						values.push_back(static_cast<integer>(nodes.size()));
						for (ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
							values.push_back(jt->IntegerDF(set_id));
						break;
					}
					case Element::Tri:
					case Element::Quad:
					case Element::MultiLine:
					case Element::Polygon:
					{
						ElementArray<Node> nodes = it->getNodes();
						values.push_back(static_cast<integer>(nodes.size()));
						for (ElementArray<Node>::iterator jt = nodes.begin(); jt != nodes.end(); jt++)
							values.push_back(jt->IntegerDF(set_id));
						break;
					}
					default: std::cout << __FILE__ << ":" << __LINE__ << " This should not happen " << Element::GeometricTypeName(it->GetGeometricType()) << std::endl;
					}
				}
				f << "CELLS " << NumberOfCells() + NumberOfFaces() << " " << values.size() << "\n";
			}
			else
				f << "CELLS " << NumberOfCells() << " " << values.size() << "\n";

			for(size_t i = 0; i < values.size(); i++)
			{
				f << values[i] << " ";
				if( (i+1) % 20 == 0) f << "\n";
			}
			f << "\n";
		}
		if (output_faces)
			f << "CELL_TYPES " << NumberOfCells() + NumberOfFaces() << "\n";
		else 			
			f << "CELL_TYPES " << NumberOfCells() << "\n";

		for(Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
		{
			INMOST_DATA_ENUM_TYPE nnodes = VtkElementNodes(it->GetGeometricType());
			if( nnodes == ENUMUNDEF || nnodes == it->nbAdjElements(NODE) ) //nodes match - output correct type
				f << VtkElementType(it->GetGeometricType()) << "\n";
			else //number of nodes mismatch with expected - some topology checks must be off
				f << VtkElementType(Element::MultiPolygon) << "\n";
		}
		if (output_faces){
			for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); it++)
			{
				INMOST_DATA_ENUM_TYPE nnodes = VtkElementNodes(it->GetGeometricType());
				if (nnodes == ENUMUNDEF || nnodes == it->nbAdjElements(NODE)) //nodes match - output correct type
					f << VtkElementType(it->GetGeometricType()) << "\n";
				else //number of nodes mismatch with expected - some topology checks must be off
					f << VtkElementType(Element::MultiPolygon) << "\n";
			}
		}
		DeleteTag(set_id);
		{
			std::vector<std::string> tag_names;
			std::vector<Tag> tags;
			ListTagNames(tag_names);
			for(unsigned int i = 0; i < tag_names.size(); i++)
			{
				Tag t = GetTag(tag_names[i]);
				if (((t.isDefined(CELL) /*&& !t.isSparse(CELL)*/)
					|| (t.isDefined(FACE) && output_faces)) &&
						t.GetPrint() && //Temporary solution: @see Mesh::file_option
						t.GetDataType() != DATA_BULK && 
						t.GetDataType() != DATA_REFERENCE &&
						t.GetDataType() != DATA_REMOTE_REFERENCE &&
						t != CoordsTag() && 
						t != SharedTag() && 
						t != SendtoTag() && 
						t != ProcessorsTag())
				{
					bool skip = CheckSaveSkip(tag_names[i],nosave,saveonly);
					if( !skip )
						tags.push_back(t);
				}
			}
				
			if (!tags.empty() && output_faces) 
			{ 
				f << "CELL_DATA " << NumberOfCells() + NumberOfFaces() << "\n"; 
			}
			else if (!tags.empty()) f << "CELL_DATA " << NumberOfCells() << "\n";

			for(INMOST_DATA_ENUM_TYPE i = 0; i < tags.size(); i++)
			{
				INMOST_DATA_ENUM_TYPE comps = tags[i].GetSize();
				if( comps == ENUMUNDEF )
					continue;
				else
				{
					{
						std::string type_str = "int";
						if(  tags[i].GetDataType() == DATA_REAL
#if defined(USE_AUTODIFF)
						   || tags[i].GetDataType() == DATA_VARIABLE
#endif
						   ) type_str = "double";
						f << "SCALARS " << Space2Underscore(tags[i].GetTagName()) << " " << type_str << " " << comps << "\n";
						f << "LOOKUP_TABLE default\n";
						for (Mesh::iteratorCell it = BeginCell(); it != EndCell(); it++)
						{
							switch (tags[i].GetDataType())
							{
							case DATA_REAL:
							{
								if (tags[i].isDefined(CELL) && it->HaveData(tags[i]))
								{
									Storage::real_array arr = it->RealArray(tags[i]);
									for (unsigned int m = 0; m < comps; m++)
									{
										double val = static_cast<double>(arr[m]);
										f << (__isbad(val) ? -0.9999E30 : val) << " ";
									}
								}
								else for (unsigned int m = 0; m < comps; m++) f << -0.9999E30 << " ";
								f << "\n";
							}
							break;
							case DATA_INTEGER:
							{
								if (tags[i].isDefined(CELL) && it->HaveData(tags[i]))
								{
									Storage::integer_array arr = it->IntegerArray(tags[i]);
									for (unsigned int m = 0; m < comps; m++) f << arr[m] << " ";
								}
								else for (unsigned int m = 0; m < comps; m++) f << INT_MIN << " ";
								f << "\n";
							}
							break;
#if defined(USE_AUTODIFF)
							case DATA_VARIABLE:
							{
								if (tags[i].isDefined(CELL) && it->HaveData(tags[i]))
								{
									Storage::var_array arr = it->VariableArray(tags[i]);
									for (unsigned int m = 0; m < comps; m++)
									{
										double val = static_cast<double>(arr[m].GetValue());
										f << (__isbad(val) ? -0.9999E30 : val) << " ";
									}
								}
								else for (unsigned int m = 0; m < comps; m++) f << -0.9999E30 << " ";
								f << "\n";
							}
							break;
#endif
							default: continue;
							}
						}
							


						if (output_faces)
						{
							for (Mesh::iteratorFace it = BeginFace(); it != EndFace(); it++)
							{
								switch (tags[i].GetDataType())
								{
								case DATA_REAL:
								{
									if (tags[i].isDefined(FACE) && it->HaveData(tags[i]))
									{
										Storage::real_array arr = it->RealArray(tags[i]);
										for (unsigned int m = 0; m < comps; m++) f << arr[m] << " ";
									}
									else for (unsigned int m = 0; m < comps; m++) f << -0.9999E30 << " ";
									f << "\n";
								}
								break;
								case DATA_INTEGER:
								{
									if (tags[i].isDefined(FACE) && it->HaveData(tags[i]))
									{
										Storage::integer_array arr = it->IntegerArray(tags[i]);
										for (unsigned int m = 0; m < comps; m++) f << arr[m] << " ";
									}
									else for (unsigned int m = 0; m < comps; m++) f << INT_MIN << " ";
									f << "\n";
								}
								break;
#if defined(USE_AUTODIFF)
								case DATA_VARIABLE:
								{
									if (tags[i].isDefined(FACE) && it->HaveData(tags[i]))
									{
										Storage::var_array arr = it->VariableArray(tags[i]);
										for (unsigned int m = 0; m < comps; m++)
										{
											double val = static_cast<double>(arr[m].GetValue());
											f << (__isbad(val) ? -0.9999E30 : val) << " ";
										}
									}
									else for (unsigned int m = 0; m < comps; m++) f << -0.9999E30 << " ";
									f << "\n";
								}
								break;
#endif
								default: continue;
								}
							}
						}
					}
				}
			}
		}
		{
			std::vector<std::string> tag_names;
			std::vector<Tag> tags;
			ListTagNames(tag_names);
			for(unsigned int i = 0; i < tag_names.size(); i++)
			{
				Tag t = GetTag(tag_names[i]);
				if( t.isDefined(NODE) && /*!t.isSparse(NODE) &&*/
				    t.GetDataType() != DATA_BULK &&
				    t.GetDataType() != DATA_REFERENCE &&
				    t.GetDataType() != DATA_REMOTE_REFERENCE &&
					t != CoordsTag() &&
				    t != SharedTag() &&
				    t != SendtoTag() &&
				    t != ProcessorsTag())
					tags.push_back(t);
			}
				
			if( !tags.empty() ) f << "POINT_DATA " << NumberOfNodes() << "\n";
			for(INMOST_DATA_ENUM_TYPE i = 0; i < tags.size(); i++)
			{
				INMOST_DATA_ENUM_TYPE comps = tags[i].GetSize();
				if( comps == ENUMUNDEF )
					continue;
				else
				{
					{
						std::string type_str = "int";
						if(  tags[i].GetDataType() == DATA_REAL
#if defined(USE_AUTODIFF)
						   || tags[i].GetDataType() == DATA_VARIABLE
#endif
						   ) type_str = "double";
						f << "SCALARS " << Space2Underscore(tags[i].GetTagName()) << " " << type_str << " " << comps << "\n";
						f << "LOOKUP_TABLE default\n";
						for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); it++)
						{
							switch( tags[i].GetDataType() )
							{
								case DATA_REAL:
								{
									if (it->HaveData(tags[i]))
									{
										Storage::real_array arr = it->RealArray(tags[i]);
										for(unsigned int m = 0; m < comps; m++)
										{
											double val = static_cast<double>(arr[m]);
											f << (__isbad(val) ? -0.9999E30 : val) << " ";
										}
										f << "\n";
									}
									else for(unsigned int m = 0; m < comps; m++) f << -0.9999E30 << " ";
								}
								break;
								case DATA_INTEGER:
								{
									if (it->HaveData(tags[i]))
									{
										Storage::integer_array arr = it->IntegerArray(tags[i]);
										for(unsigned int m = 0; m < comps; m++) f << arr[m] << " ";
										f << "\n";
									}
									else for (unsigned int m = 0; m < comps; m++) f << INT_MIN << " ";
								}
								break;
#if defined(USE_AUTODIFF)
								case DATA_VARIABLE:
								{
									if (it->HaveData(tags[i]))
									{
										Storage::var_array arr = it->VariableArray(tags[i]);
										for(unsigned int m = 0; m < comps; m++)
										{
											double val = static_cast<double>(arr[m].GetValue());
											f << (__isbad(val) ? -0.9999E30 : val) << " ";
										}
										f << "\n";
									}
									else for(unsigned int m = 0; m < comps; m++) f << -0.9999E30 << " ";
								}
								break;
#endif
								default: continue;
							}
						}
					}
				}
			}
		}
	}

	void Mesh::LoadVTK(std::string File)
	{
		int verbosity = 0;
		for(INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
		{
			if( file_options[k].first == "VERBOSITY" )
			{
				verbosity = atoi(file_options[k].second.c_str());
				if( verbosity < 0 || verbosity > 2 )
				{
					std::cout << __FILE__ << ":" << __LINE__ << " Unknown verbosity option: " << file_options[k].second << std::endl;
					verbosity = 1;
				}
			}
		}
		
		std::set< std::string > noload, loadonly;
		
		noload = TagOptions("noload");
		loadonly = TagOptions("loadonly");
		

		MarkerType unused_marker = CreateMarker();
		int grid_is_2d = 2;
		for(INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); ++k)
		{
			if( file_options[k].first == "VTK_GRID_DIMS" )
			{
				if( file_options[k].second == "AUTO" )
					grid_is_2d = 2;
				if( atoi(file_options[k].second.c_str()) == 2 )
					grid_is_2d = 1;
				else if( atoi(file_options[k].second.c_str()) == 3 )
					grid_is_2d = 0;
			}
		}

		//Determine whether there are already 3d elements so that the grid is 3d
		if( grid_is_2d == 2 && NumberOfCells() )
		{
			for(Mesh::iteratorCell it = BeginCell(); it != EndCell() && grid_is_2d == 2; ++it)
				if( it->GetElementDimension() == 3 )
					grid_is_2d = 0;
		}


		std::vector<HandleType> old_nodes(NumberOfNodes());
		{
			unsigned qq = 0;
			for(Mesh::iteratorNode it = BeginNode(); it != EndNode(); ++it)
				old_nodes[qq++] = *it;
		}
		if( !old_nodes.empty() ) 
		{
			std::sort(old_nodes.begin(),old_nodes.end(),CentroidComparator(this));
			//for(std::vector<HandleType>::iterator it = old_nodes.begin(); it != old_nodes.end(); ++it)
			//{
			//	Storage::real_array c = RealArrayDF(*it,CoordsTag());
			//	REPORT_VAL("coord: ",c[0] << " " << c[1] << " " << c[2]);
			//}
		}

			
		FILE * f = fopen(File.c_str(),"r");
		if( !f ) 
		{
			std::cout << __FILE__ << ":" << __LINE__ << " cannot open " << File << std::endl;
			throw BadFileName;
		}
		std::vector<Tag> datatags;
		std::vector<HandleType> newnodes;
		std::vector<HandleType> newcells;
		std::vector<int> cp;
		std::vector<int> ct;
		unsigned int state = R_VERSION;
		char readline[2048];
		int filled;
		bool binary = false;
		int read_into = NONE, read_into_cell = NONE;
		bool read_next_line = true;
		while( state != R_QUIT )
		{
			
			if( read_next_line )
			{
				if( fgets(readline,2048,f) == NULL )
				{
					state = R_QUIT;
					continue;
				}
			} else read_next_line = true;
			//if( readline[strlen(readline)-1] == '\n' )
			//	readline[strlen(readline)-1] = '\0';
			while( strlen(readline) != 0 && isspace(readline[strlen(readline)-1]) )
				readline[strlen(readline)-1] = '\0';
			if( strlen(readline) == 0 ) continue;
			switch( state )
			{
				case R_VERSION:
				{
					int h,l;
					filled = sscanf(readline,"# vtk DataFile Version %1d.%1d",&h,&l);
					if( filled == 0 ) 
					{
						std::cout << __FILE__ << ":" << __LINE__ << " version information not found in " << File << std::endl;
						throw BadFile;
					}
					//file version checker
					state = R_USERDATA;
					break;
				}
				case R_USERDATA:
				{
					if (!strncmp(readline, "VTK_OUTPUT_FACES", 16)) this->SetFileOption("VTK_OUTPUT_FACES", "1");
					state = R_DATATYPE;
					break;
				}
				case R_DATATYPE:
				{
					if( !strncmp(readline,"BINARY",6) )
					{
						binary = 1;
						state = R_WAITDATA;
					}
					else if( !strncmp(readline,"ASCII",5) )
						state = R_WAITDATA;
					else 
					{
						std::cout << __FILE__ << ":" << __LINE__ << " unexpected data type " << readline << " in " << File << " expected BINARY or ASCII instead " << std::endl;
						throw BadFile;
					}
					break;
				}
				case R_ATTRDATA:
				{
					char dataname[1024];
					char attrname[1024];
					char attrtype[1024];
					int nentries = 1;
					filled = sscanf(readline," %s ",dataname);
					if( filled != 1 ) 
					{
						std::cout << __FILE__ << ":" << __LINE__ << " cannot read attribute name in " << File << std::endl;
						throw BadFile;
					}
					if( !strcmp(dataname,"SCALARS") )
					{
						DataType t = DATA_BULK;
						filled = sscanf(readline," SCALARS %s %s %d",attrname,attrtype,&nentries);
						bool skip = CheckLoadSkip(attrname,noload,loadonly);
						if( verbosity > 0 ) std::cout << (skip?"Skipping":"Reading") << " scalar attribute " << attrname << std::endl;
						if( filled < 2 ) 
						{
							std::cout << __FILE__ << ":" << __LINE__ << " found " << filled << " arguments to SCALARS field, must be >= 2\nline:\n" << readline << "\n";
							throw BadFile;
						}
						for(unsigned int i = 0; i < strlen(attrtype); i++) attrtype[i] = tolower(attrtype[i]);
						if( !strcmp(attrtype,"bit") || !strcmp(attrtype,"unsigned_char") || !strcmp(attrtype,"char") ||
							  !strcmp(attrtype,"unsigned_short") || !strcmp(attrtype,"short") || !strcmp(attrtype,"unsigned_int") ||
							  !strcmp(attrtype,"int") || !strcmp(attrtype,"unsigned_long") || !strcmp(attrtype,"long"))
							t = DATA_INTEGER;
						else if( !strcmp(attrtype,"float") || !strcmp(attrtype,"double") )
							t = DATA_REAL;
						else
						{
							std::cout << __FILE__ << ":" << __LINE__ << " unexpected data type " << attrtype << " in " << File << " for " << attrname << std::endl;
							std::cout << "expected one of: bit, unsigned_char, char, unsigned_short, short, unsigned_int, int, unsigned_long, long, float, double" << std::endl;
							throw BadFile;
						}
						if( fgets(readline,2048,f) == NULL ) 
						{
							std::cout << __FILE__ << ":" << __LINE__ << " expected LOOKUP_TABLE in " << File << std::endl;
							throw BadFile; //LOOK_UP TABLE
						}
						if( read_into == 2 )
						{
							bool sparse = (std::string(attrname).substr(0,9) != "PROTECTED");
							unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
							if( skip )
							{
								for(unsigned int it = 0; it < newcells.size(); it++)
								{

									if( t == DATA_INTEGER )
										for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile;}
									if( t == DATA_REAL )
										for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; }
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newcells.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							else
							{
								Tag attr = CreateTag(attrname,t,read_into_cell,sparse ? read_into_cell & ~CELL : NONE,nentries);
								for(unsigned int it = 0; it < newcells.size(); it++)
								{

									if( t == DATA_INTEGER )
									{
										if( newcells[it] != InvalidHandle() )
										{
											Storage::integer_array attrdata = IntegerArray(newcells[it],attr);
											for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
										}
										else for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile;}
									}
									if( t == DATA_REAL )
									{
										if( newcells[it] != InvalidHandle() )
										{
											Storage::real_array attrdata = RealArray(newcells[it],attr);
											for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
										}
										else for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; }
									}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newcells.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
								datatags.push_back(attr);
							}
						}
						if( read_into == 1 )
						{
							unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size())/250,1);
							if( skip )
							{
								for(unsigned int it = 0; it < newnodes.size(); it++)
								{
									if( t == DATA_INTEGER )
										for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile;}
									if( t == DATA_REAL )
										for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile;}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newnodes.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							else
							{
								Tag attr = CreateTag(attrname,t,NODE,NONE,nentries);
								for(unsigned int it = 0; it < newnodes.size(); it++)
								{
									if( t == DATA_INTEGER )
									{
										Storage::integer_array attrdata = IntegerArray(newnodes[it],attr);
										for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
									}
									if( t == DATA_REAL )
									{
										Storage::real_array attrdata = RealArray(newnodes[it],attr);
										for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
									}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newnodes.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							//filled = fscanf(f,"\n");
						}
						if( fgets(readline,2048,f) == NULL ) throw BadFile;
						skip_empty(readline,f);
						skip_metadata(3,readline,f);
						skip_empty(readline,f);
						read_next_line = false;
						
						break;
					}
					else if( !strcmp(dataname,"COLOR_SCALARS") )
					{	
						char attrname[1024];
						filled = sscanf(readline," COLOR_SCALARS %s %d",attrname,&nentries);
						bool skip = CheckLoadSkip(attrname,noload,loadonly);
						if( verbosity > 0 ) std::cout << (skip?"Skipping":"Reading") << " color attribute " << attrname << "\n";
						if( read_into == 2 )
						{
							bool sparse = (std::string(attrname).substr(0,9) != "PROTECTED");
							unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
							if( skip )
							{
								for(unsigned int it = 0; it < newcells.size(); it++)
								{
									for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile;}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newcells.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							else
							{
								Tag attr = CreateTag(attrname,DATA_REAL,read_into_cell,sparse ? read_into_cell & ~CELL : NONE,nentries);
								for(unsigned int it = 0; it < newcells.size(); it++)
								{
									if( newcells[it] != InvalidHandle() )
									{
										Storage::real_array attrdata = RealArray(newcells[it],attr);
										for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
									}
									else for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile;}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newcells.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
								datatags.push_back(attr);
							}
						}
						if( read_into == 1 )
						{
							unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size()/250),1);
							if( skip ) 
							{
								for(unsigned int it = 0; it < newnodes.size(); it++)
								{
									for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile;}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newnodes.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							else
							{
								Tag attr = CreateTag(attrname,DATA_REAL,NODE,NONE,nentries);
								for(unsigned int it = 0; it < newnodes.size(); it++)
								{
									Storage::real_array attrdata = RealArray(newnodes[it],attr);
									for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newnodes.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							//filled = fscanf(f,"\n");
						}
						if( fgets(readline,2048,f) == NULL ) throw BadFile;
						skip_empty(readline,f);
						skip_metadata(3,readline,f);
						skip_empty(readline,f);
						read_next_line = false;
					
						break;
					}
					else if( !strcmp(dataname," LOOKUP_TABLE") )
					{
						int nentries;
						filled = sscanf(readline," LOOKUP_TABLE %*s %d",&nentries);
						if( filled != 1 ) throw BadFile;
						for(int i = 0; i < nentries; i++)
						{
							filled = fscanf(f," %*f %*f %*f %*f");
						}
						filled = fscanf(f,"\n");
						break;
					}
					else if( !strcmp(dataname,"NORMALS") || !strcmp(dataname,"VECTORS") || !strcmp(dataname,"TENSORS"))
					{
						char attrname[1024];
						char attrtype[1024];
						int nentries = 3;
						if( !strcmp(dataname,"TENSORS") ) nentries = 9;
						DataType t = DATA_BULK;
						filled = sscanf(readline,"%*s %s %s",attrname,attrtype);
						bool skip = CheckLoadSkip(attrname,noload,loadonly);
						if( verbosity > 0 ) std::cout << (skip? "Skipping":"Reading") << " " << dataname << " attribute " << attrname << "n";
						if( filled != 2 ) throw BadFile;
						for(unsigned int i = 0; i < strlen(attrtype); i++) attrtype[i] = tolower(attrtype[i]);
						if( !strcmp(attrtype,"bit") || !strcmp(attrtype,"unsigned_char") || !strcmp(attrtype,"char") ||
							  !strcmp(attrtype,"unsigned_short") || !strcmp(attrtype,"short") || !strcmp(attrtype,"unsigned_int") ||
							  !strcmp(attrtype,"int") || !strcmp(attrtype,"unsigned_long") || !strcmp(attrtype,"long"))
							t = DATA_INTEGER;
						else if( !strcmp(attrtype,"float") || !strcmp(attrtype,"double") )
							t = DATA_REAL;
						else
							throw BadFile;
						if( fgets(readline,2048,f) == NULL ) throw BadFile; //LOOK_UP TABLE
						if( read_into == 2 )
						{
							bool sparse = (std::string(attrname).substr(0,9) != "PROTECTED");
							unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
							if( skip )
							{
								for(unsigned int it = 0; it < newcells.size(); it++)
								{
									if( t == DATA_INTEGER )
										for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile;}
									if( t == DATA_REAL )
										for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile;}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newcells.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
							}
							else
							{
								Tag attr = CreateTag(attrname,t,read_into_cell,sparse ? read_into_cell & ~CELL : NONE,nentries);
								for(unsigned int it = 0; it < newcells.size(); it++)
								{
									if( t == DATA_INTEGER )
									{
										if( newcells[it] != InvalidHandle() )
										{
											Storage::integer_array attrdata = IntegerArray(newcells[it],attr);
											for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
										}
										else for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile;}
									}
									if( t == DATA_REAL )
									{
										if( newcells[it] != InvalidHandle() )
										{
											Storage::real_array attrdata = RealArray(newcells[it],attr);
											for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
										}
										else for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile;}
									}
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newcells.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}
								}
								datatags.push_back(attr);
							}
						}
						if( read_into == 1 )
						{
							unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size()/250),1);
							if( skip )
							{
								for(unsigned int it = 0; it < newnodes.size(); it++)
								{
									if( t == DATA_INTEGER )
										for(int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile; }
									if( t == DATA_REAL )
										for(int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; }
									if( verbosity > 1 && it%report_pace == 0)
									{
										std::ios save(NULL);
										save.copyfmt(std::cout);
										std::cout << "data " << std::setw(5) << std::fixed << std::setprecision(1) << (it*100.0)/(1.0*newnodes.size()) << "%\r" << std::flush;
										std::cout.copyfmt(save);
									}