mesh_file.cpp

#ifdef _MSC_VER //kill some warnings
#define _CRT_SECURE_NO_WARNINGS
#endif
#include "inmost.h"

#if defined(USE_MESH)
#include <deque>
#include <fstream>
#include <sstream>
#include <stdlib.h>
#include <stdio.h>
#include "io.hpp"


//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


//gmsh states

#define GMSH_VERNONE 0
#define GMSH_VER1 1
#define GMSH_VER2 2

#define GMSH_NONE 0
#define GMSH_NODES_TOT 1
#define GMSH_NODES_NUM 2
#define GMSH_NODES_X 3
#define GMSH_NODES_Y 4
#define GMSH_NODES_Z 5
#define GMSH_ELEMENTS_TOT 6
#define GMSH_ELEMENTS_NUM 7
#define GMSH_ELEMENTS_TYPE 8
#define GMSH_ELEMENTS_REGPHYS 9
#define GMSH_ELEMENTS_REGELEM 10
#define GMSH_ELEMENTS_NUMTAGS 11
#define GMSH_ELEMENTS_TAGS 12
#define GMSH_ELEMENTS_NUMNODES 13
#define GMSH_ELEMENTS_NODELIST 14
#define GMSH_FORMAT 15
#define GMSH_SKIP_KEYWORD 16

//eclipse states
#define ECLSTRCMP(x,y) strncmp(x,y,8)

#define ECL_NEVER -1
#define ECL_NONE 0
#define ECL_SKIP_SECTION 1
#define ECL_INCLUDE 2
#define ECL_DIMENS 3
#define ECL_DX 4
#define ECL_DY 5
#define ECL_DZ 6
#define ECL_TOPS 7
#define ECL_PERMX 8
#define ECL_PERMY 9
#define ECL_PERMZ 10
#define ECL_PORO 11
#define ECL_MAPAXIS 12
#define ECL_INRAD 13
#define ECL_COORDS 14
#define ECL_ZCORN 15

#define ECL_GTYPE_NONE 0
#define ECL_GTYPE_TOPS 1
#define ECL_GTYPE_ZCORN 2
#define ECL_GTYPE_RADIAL 3
#define ECL_GTYPE_CARTESIAN 4


#define ECL_VAR_NONE 0
#define ECL_VAR_REAL 1
#define ECL_VAR_INT 2

#define ECL_IJK_DATA(i,j,k) (i + (j+k*dims[1])*dims[0])



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;
	}
}



#if defined(USE_PARALLEL_WRITE_TIME)
#define REPORT_MPI(x) {WriteTab(out_time) << "<MPI><![CDATA[" << #x << "]]></MPI>\n"; x;}
#define REPORT_STR(x) {WriteTab(out_time) << "<TEXT><![CDATA[" << x << "]]></TEXT>\n";}
#define REPORT_VAL(str,x) {WriteTab(out_time) << "<VALUE name=\"" << str << "\"> <CONTENT><![CDATA[" << x << "]]> </CONTENT><CODE><![CDATA[" << #x << "]]></CODE></VALUE>\n";}
#define ENTER_FUNC() long double all_time = Timer(); WriteTab(out_time) << "<FUNCTION name=\"" << __FUNCTION__ << "\" id=\"func" << func_id++ << "\">\n"; Enter();
#define EXIT_FUNC() WriteTab(out_time) << "<TIME>" << Timer() - all_time << "</TIME>\n"; Exit(); WriteTab(out_time) << "</FUNCTION>\n";
#else
#define REPORT_MPI(x) x
#define REPORT_STR(x) 
#define REPORT_VAL(str,x)
#define ENTER_FUNC()
#define EXIT_FUNC() 
#endif


#define GMV_HEADER 0
#define GMV_NODES 1
#define GMV

namespace INMOST
{

	typedef char HeaderType;
	const HeaderType EndOfData  = 0x01;
	const HeaderType NodeHeader = 0x02;
	const HeaderType EdgeHeader = 0x03;
	const HeaderType FaceHeader = 0x04;
	const HeaderType CellHeader = 0x05;
	const HeaderType ESetHeader = 0x06;
	const HeaderType TagsHeader = 0x07;
	const HeaderType MeshHeader = 0x08;
	const HeaderType EoMHeader  = 0x09;
	const HeaderType INMOSTFile   = 0x10;
	const HeaderType MeshDataHeader = 0x11;

	void Mesh::SetFileOption(std::string key, std::string val)
	{
		for(INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); k++)
		{
			if(file_options[k].first == key)
			{
				file_options[k].second = val;
			}
		}
		file_options.push_back(std::make_pair(key,val));
	}

	std::string Mesh::GetFileOption(std::string key)
	{
		for(INMOST_DATA_ENUM_TYPE k = 0; k < file_options.size(); k++)
		{
			if(file_options[k].first == key)
			{
				return file_options[k].second;
			}
		}
		return "";
	}
	
	std::ostream & operator <<(std::ostream & out, HeaderType H)
	{
		out.put(H);
		return out;
	}
	
	std::istream & operator >>(std::istream & in, HeaderType &H)
	{
		in.get(H);
		return in;
	}

	
	void Mesh::Load(std::string File)
	{
		ENTER_FUNC();
		REPORT_VAL("File",File);
		std::string LFile;
		LFile.resize(File.size());
		std::transform(File.begin(),File.end(),LFile.begin(),::tolower);
		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 )
				{
					printf("%s:%d Unknown verbosity option: %s\n",__FILE__,__LINE__,file_options[k].second.c_str());
					verbosity = 1;
				}
			}
		}
//		long double load_timer = Timer();
		/*
		if(File.find("gmv") != std::string::npos) //this is gmv
		{
			char readline[2048];
			std::vector<File *> fs;
			fs.push_back(fopen(File.c_str(),"r"));
			if( fs[0] == NULL ) throw BadFileName;
			while( !fs.empty() )
			{
				while( fgets(readline,2048,fs[fs.size()-1]) != NULL )
				{
					switch(state)
					{
						case 
					}
				}
				fclose(fs[fs.size()-1]);
				fs.pop();
			}
		}
		else 
		*/
		if(LFile.find(".grdecl") != std::string::npos) // this is eclipse grid
		{
			FILE * f = fopen(LFile.c_str(),"r");
			if( f == NULL )
			{
				std::cout << __FILE__ << ":" << __LINE__ << " cannot open file " << LFile << std::endl;
				throw BadFileName;
			}
			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));
			std::vector< std::pair< std::pair<FILE *,std::string>, int> > fs(1,std::make_pair(std::make_pair(f,File),0));
			char readline[2048], *p, *pend, rec[2048];
			int text_end, text_start, state = ECL_NONE, nchars;
			int waitlines = 0;
			int have_dimens = 0, totread, downread, numrecs, offset;
			int gtype = ECL_GTYPE_NONE;
			int argtype = ECL_VAR_NONE;
			int radial = ECL_GTYPE_NONE;
			Storage::real * read_arrayf = NULL;
			Storage::integer * read_arrayi = NULL;
			Storage::integer dims[3], mapaxis[6] = {0,1,0,0,1,0};
			Storage::real inrad = 0;
			std::vector<Storage::real> xyz,perm,poro, tops,zcorn;
			std::vector<Storage::integer> actnum;
			while(!fs.empty())
			{
				while(fgets(readline,2048,fs.back().first.first) != NULL)
				{
					fs.back().second++; //line number
					{
						if( readline[strlen(readline)-1] == '\n' ) readline[strlen(readline)-1] = '\0';
						text_end = static_cast<int>(strlen(readline));
						for(text_start = 0; isspace(readline[text_start]) && text_start < text_end; text_start++);
						if( text_start == text_end ) continue;
						for(text_end = text_end-1; isspace(readline[text_end]) && text_end > text_start; text_end--);
						readline[text_end+1] = '\0';
						p = readline + text_start;
						pend = readline + text_end + 1;
						for(char * q = p; q < pend; q++) *q = toupper(*q);
					}
					if( p[0] == '-' && p[1] == '-' ) continue; //skip comment
					if(waitlines) {waitlines--; continue;} //skip meaningful lines
					switch(state)
					{
					case ECL_NONE:
						if( !ECLSTRCMP(p,"END") ) //end of data - don't read beyond
						{
							goto ecl_exit_loop;
						}
						else if( !ECLSTRCMP(p,"INCLUDE") ) state = ECL_INCLUDE;
						else if( !ECLSTRCMP(p,"DIMENS") || !ECLSTRCMP(p,"SPECGRID") )
						{
							read_arrayi = dims;
							numrecs = 1;
							downread = totread = 3;
							argtype = ECL_VAR_INT;
							offset = state = ECL_DIMENS;
							have_dimens = 1;
						}
						else if( !ECLSTRCMP(p,"MAPAXIS") )
						{
							read_arrayi = mapaxis;
							numrecs = 1;
							downread = totread = 6;
							argtype = ECL_VAR_INT;
							offset = state = ECL_MAPAXIS;
						}
						else if( !ECLSTRCMP(p,"DX") )
						{
							assert(have_dimens);
							if( xyz.empty() ) xyz.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = xyz.empty()? NULL : &xyz[0];
							numrecs = 3;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = state = ECL_DX;
						}
						else if( !ECLSTRCMP(p,"DY") )
						{
							assert(have_dimens);
							if( xyz.empty() ) xyz.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = xyz.empty() ? NULL : &xyz[0];
							numrecs = 3;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = ECL_DX;
							state = ECL_DY;
						}
						else if( !ECLSTRCMP(p,"DZ") )
						{
							assert(have_dimens);
							if( xyz.empty() ) xyz.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = xyz.empty() ? NULL : &xyz[0];
							numrecs = 3;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = ECL_DX;
							state = ECL_DZ;
						}
						else if( !ECLSTRCMP(p,"COORD") )
						{
							assert(have_dimens);
							if( xyz.empty() ) xyz.resize(2*(dims[0]+1)*(dims[1]+1));
							read_arrayf = xyz.empty() ? NULL : &xyz[0];
							numrecs = 1;
							downread = totread = 2*(dims[0]+1)*(dims[1]+1);
							argtype = ECL_VAR_REAL;
							offset = state = ECL_COORDS;
							gtype = ECL_GTYPE_ZCORN;
						}
						else if( !ECLSTRCMP(p,"ZCORN") )
						{
							assert(have_dimens);
							if( zcorn.empty() ) zcorn.resize(dims[0]*dims[1]*dims[2]*8);
							read_arrayf = zcorn.empty() ? NULL : &zcorn[0];
							numrecs = 1;
							downread = totread = dims[0]*dims[1]*dims[2]*8;
							argtype = ECL_VAR_REAL;
							state = offset = ECL_ZCORN;
						}
						else if( !ECLSTRCMP(p,"TOPS") )
						{
							assert(have_dimens);
							tops.resize(dims[0]*dims[1]*dims[2]);
							read_arrayf = tops.empty() ? NULL : &tops[0];
							numrecs = 1;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = state = ECL_TOPS;
							gtype = ECL_GTYPE_TOPS;
						}
						else if( !ECLSTRCMP(p,"PERMX") )
						{
							assert(have_dimens);
							if( perm.empty() ) perm.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = perm.empty() ? NULL : &perm[0];
							numrecs = 3;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = state = ECL_PERMX;
						}
						else if( !ECLSTRCMP(p,"PERMY") )
						{
							assert(have_dimens);
							if( perm.empty() ) perm.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = perm.empty() ? NULL : &perm[0];
							numrecs = 3;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = ECL_PERMX;
							state = ECL_PERMY;
						}
						else if( !ECLSTRCMP(p,"PERMZ") )
						{
							assert(have_dimens);
							if( perm.empty() ) perm.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = perm.empty() ? NULL : &perm[0];
							numrecs = 3;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = ECL_PERMX;
							state = ECL_PERMZ;
						}
						else if( !ECLSTRCMP(p,"PORO") )
						{
							assert(have_dimens);
							poro.resize(3*dims[0]*dims[1]*dims[2]);
							read_arrayf = poro.empty() ? NULL : &poro[0];
							numrecs = 1;
							downread = totread = dims[0]*dims[1]*dims[2];
							argtype = ECL_VAR_REAL;
							offset = state = ECL_PORO;
						}
						else if( !ECLSTRCMP(p,"RADIAL") )
						{
							radial = ECL_GTYPE_RADIAL;
						}
						else if( !ECLSTRCMP(p,"CART") )
						{
							radial = ECL_GTYPE_CARTESIAN;
						}
						else if( !ECLSTRCMP(p,"INRAD") )
						{
							if( radial != ECL_GTYPE_RADIAL ) 
							{
								std::cout << __FILE__ << ":" << __LINE__ << " inner radius specified for cartesian grid ";
								std::cout << " in " << fs.back().first.second << ":" << fs.back().second << std::endl;
							}
							if( radial == ECL_GTYPE_NONE ) radial = ECL_GTYPE_RADIAL;
							state = ECL_INRAD;
						}
						else
						{
							std::cout << __FILE__ << ":" << __LINE__ << " skipped " << p << " in " << fs.back().first.second << ":" << fs.back().second << std::endl;
						}
						break;
					case ECL_SKIP_SECTION:
						if( *(pend-1) == '/' || *p == '/' ) state = ECL_NONE;
						break;
					case ECL_INCLUDE:
						if( 1 == sscanf(p," %s",rec) )
						{
							int shift_one = 0;
							if( (rec[0] == '\'' || rec[0] == '"') && rec[0] == rec[strlen(rec)-1] ) //remove quotes
							{
								rec[strlen(rec)-1] = '\0';
								shift_one = 1;
							}
							f = fopen(rec+shift_one,"r");
							if( f == NULL )
							{
								std::cout << __FILE__ << ":" << __LINE__ << " cannot open file " << rec+shift_one << " included from "<< fs.back().first.second << " line " << fs.back().second << std::endl;
								throw BadFileName;
							}
							fs.push_back(std::make_pair(std::make_pair(f,std::string(rec+shift_one)),0));
							if( *(pend-1) == '/' ) state = ECL_NONE; else state = ECL_SKIP_SECTION;
						}
						else
						{
							std::cout << __FILE__ << ":" << __LINE__ << " cannot read file name, string " << p << " in " << fs.back().first.second << " line " << fs.back().second << std::endl;
							throw BadFile;
						}
						break;
					case ECL_ZCORN:
					case ECL_COORDS:
					case ECL_MAPAXIS:
					case ECL_DIMENS:
					case ECL_DX:
					case ECL_DY:
					case ECL_DZ:
					case ECL_PERMX:
					case ECL_PERMY:
					case ECL_PERMZ:
					case ECL_PORO:
					case ECL_TOPS:
						while( downread > 0 && p < pend )
						{
							if( 1 == sscanf(p,"%s%n",rec,&nchars) )
							{
								p += nchars;
								while(isspace(*p) && p < pend) ++p;
								int count = 1;
								int begval = 0;
								for(int q = 0; q < static_cast<int>(strlen(rec)); ++q)
									if( rec[q] == '*' )
									{
										begval = q+1;
										rec[q] = '\0';
										break;
									}
								if( begval > 0 ) count = atoi(rec);
								if( argtype == ECL_VAR_REAL )
								{
									Storage::real val = atof(rec+begval);
									while(count)
									{
										read_arrayf[numrecs*(totread-(downread--))+(state-offset)] = val;
										count--;
									}
								}
								else if( argtype == ECL_VAR_INT )
								{
									Storage::integer val = atoi(rec+begval);
									while(count)
									{
										read_arrayi[numrecs*(totread-(downread--))+(state-offset)] = val;
										count--;
									}
								}
								else
								{
									std::cout << __FILE__ << ":" << __LINE__ << " probably forgot to set up argument type to read " << std::endl;
									throw Impossible;
								}
							}
							else
							{
								std::cout << __FILE__ << ":" << __LINE__ << " cannot read data " << p << " in " << fs.back().first.second << ":" << fs.back().second << std::endl;
								throw BadFile;
							}
							if( *p == '/' ) 
							{
								if( downread > 0 )
								{
									std::cout << __FILE__ << ":" << __LINE__ << " early data termination, read " << totread-downread << " of " << totread << " records in " << fs.back().first.second << ":" << fs.back().second << std::endl;
									throw BadFile;
								}
							}
						}
						if( *(pend-1) == '/' || *p == '/' ) state = ECL_NONE; else state = ECL_SKIP_SECTION;
						break;
					case ECL_INRAD:
						if( 1 == sscanf(p,"%lf%n",&inrad,&nchars) )
						{
							p += nchars;
							while(isspace(*p) && p < pend) ++p;
						}
						else
						{
							std::cout << __FILE__ << ":" << __LINE__ << " cannot read data " << p << " in " << fs.back().first.second << ":" << fs.back().second << std::endl;
							throw BadFile;
						}
						if( *(pend-1) == '/' || *p == '/' ) state = ECL_NONE; else state = ECL_SKIP_SECTION;
						break;
					}
				}
ecl_exit_loop:
				fclose(fs.back().first.first);
				fs.pop_back();
			}
			if( radial == ECL_GTYPE_RADIAL )
			{
				std::cout << __FILE__ << ":" << __LINE__ << " radial grids not supported yet " << std::endl;
			}
			if( gtype == ECL_GTYPE_TOPS )
			{
				std::vector<HandleType> newnodes((dims[0]+1)*(dims[1]+1)*(dims[2]+1));
				Storage::real x, y, z, node_xyz[3];
				x = 0.0;
				int numnode = 0;
				for(int i = 0; i < dims[0]+1; i++)
				{
					Storage::integer pif = std::min(dims[0]-1,i), pib = std::max(i-1,0);
					y = 0.0;
					for(int j = 0; j < dims[1]+1; j++)
					{
						Storage::integer pjf = std::min(dims[1]-1,j), pjb = std::max(j-1,0);
						z = (
							 tops[ECL_IJK_DATA(pib,pjb,0)]+
							 tops[ECL_IJK_DATA(pib,pjf,0)]+
							 tops[ECL_IJK_DATA(pif,pjb,0)]+
							 tops[ECL_IJK_DATA(pif,pjf,0)]
						    )*0.25;
						z -= (
							  xyz[3*ECL_IJK_DATA(pib,pjb,0)+2]+
							  xyz[3*ECL_IJK_DATA(pib,pjf,0)+2]+
							  xyz[3*ECL_IJK_DATA(pif,pjb,0)+2]+
							  xyz[3*ECL_IJK_DATA(pif,pjf,0)+2]
							 )*0.25;
						for(int k = 0; k < dims[2]+1; k++)
						{
							Storage::integer pkf = std::min(dims[2]-1,k), pkb = std::max(k-1,0);
							node_xyz[0] = x;
							node_xyz[1] = y;
							node_xyz[2] = z;
							int find = -1;
							if( !old_nodes.empty() )
							{
								std::vector<HandleType>::iterator it = std::lower_bound(old_nodes.begin(),old_nodes.end(),node_xyz,CentroidComparator(this));
								if( it != old_nodes.end() ) 
								{
									Storage::real_array c = RealArrayDF(*it,CoordsTag());
									if( CentroidComparator(this).Compare(node_xyz,c.data()) == 0 )
										find = static_cast<int>(it - old_nodes.begin());
								}
							}
							if( find == -1 ) newnodes[numnode++] = CreateNode(node_xyz)->GetHandle();
							else newnodes[numnode++] = old_nodes[find];
							//std::cout << i << " " << j << " " << k << " ( " << x << " , " << y << " , " << z << ") " << newnodes.back()->LocalID() << std::endl; 
							x += (
								  (
								   xyz[3*ECL_IJK_DATA(pib,pjb,pkf)+0]+
								   xyz[3*ECL_IJK_DATA(pib,pjf,pkf)+0]+
								   xyz[3*ECL_IJK_DATA(pif,pjb,pkf)+0]+
								   xyz[3*ECL_IJK_DATA(pif,pjf,pkf)+0]
								  )
								-
								  (
								   xyz[3*ECL_IJK_DATA(pib,pjb,pkb)+0]+
								   xyz[3*ECL_IJK_DATA(pib,pjf,pkb)+0]+
								   xyz[3*ECL_IJK_DATA(pif,pjb,pkb)+0]+
								   xyz[3*ECL_IJK_DATA(pif,pjf,pkb)+0]
								  )
								 )*0.25;
							y += (
								  (
							       xyz[3*ECL_IJK_DATA(pib,pjb,pkf)+1]+
								   xyz[3*ECL_IJK_DATA(pib,pjf,pkf)+1]+
								   xyz[3*ECL_IJK_DATA(pif,pjb,pkf)+1]+
								   xyz[3*ECL_IJK_DATA(pif,pjf,pkf)+1]
								  )
								-
								  (
								   xyz[3*ECL_IJK_DATA(pib,pjb,pkb)+1]+
								   xyz[3*ECL_IJK_DATA(pib,pjf,pkb)+1]+
								   xyz[3*ECL_IJK_DATA(pif,pjb,pkb)+1]+
								   xyz[3*ECL_IJK_DATA(pif,pjf,pkb)+1]
								  )
								 )*0.25;
							z += (
								  xyz[3*ECL_IJK_DATA(pib,pjb,pkb)+2]+
								  xyz[3*ECL_IJK_DATA(pib,pjf,pkb)+2]+
								  xyz[3*ECL_IJK_DATA(pif,pjb,pkb)+2]+
								  xyz[3*ECL_IJK_DATA(pif,pjf,pkb)+2]+
							      xyz[3*ECL_IJK_DATA(pib,pjb,pkf)+2]+
								  xyz[3*ECL_IJK_DATA(pib,pjf,pkf)+2]+
								  xyz[3*ECL_IJK_DATA(pif,pjb,pkf)+2]+
								  xyz[3*ECL_IJK_DATA(pif,pjf,pkf)+2]
								 )*0.125;
						}
						y += (
							  xyz[3*ECL_IJK_DATA(pib,pjb,0)+1]+
							  xyz[3*ECL_IJK_DATA(pif,pjb,0)+1]+
							  xyz[3*ECL_IJK_DATA(pib,pjf,0)+1]+
							  xyz[3*ECL_IJK_DATA(pif,pjf,0)+1]
							 )*0.25; 
					}
					x += (
						  xyz[3*ECL_IJK_DATA(pib,0,0)+0]+
						  xyz[3*ECL_IJK_DATA(pif,0,0)+0]
						 )*0.5; 
				}
				Tag tagporo,tagperm;
				if( !poro.empty() ) tagporo = CreateTag("PORO",DATA_REAL,CELL,NONE,1);
				if( !perm.empty() ) tagperm = CreateTag("PERM",DATA_REAL,CELL,NONE,3);

				const Storage::integer nvf[24] = { 2, 3, 1, 0, 4, 5, 7, 6, 0, 1, 5, 4, 3, 2, 6, 7, 2, 0, 4, 6, 1, 3, 7, 5 };
				const Storage::integer numnodes[6] = { 4, 4, 4, 4, 4, 4 };
				for(int i = 0; i < dims[0]; i++)
					for(int j = 0; j < dims[1]; j++)
						for(int k = 0; k < dims[2]; k++)
						{
							HandleType verts[8];
							verts[0] = newnodes[((i+0)*(dims[1]+1)+(j+0))*(dims[2]+1)+(k+0)];
							verts[1] = newnodes[((i+1)*(dims[1]+1)+(j+0))*(dims[2]+1)+(k+0)];
							verts[2] = newnodes[((i+0)*(dims[1]+1)+(j+1))*(dims[2]+1)+(k+0)];
							verts[3] = newnodes[((i+1)*(dims[1]+1)+(j+1))*(dims[2]+1)+(k+0)];
							verts[4] = newnodes[((i+0)*(dims[1]+1)+(j+0))*(dims[2]+1)+(k+1)];
							verts[5] = newnodes[((i+1)*(dims[1]+1)+(j+0))*(dims[2]+1)+(k+1)];
							verts[6] = newnodes[((i+0)*(dims[1]+1)+(j+1))*(dims[2]+1)+(k+1)];
							verts[7] = newnodes[((i+1)*(dims[1]+1)+(j+1))*(dims[2]+1)+(k+1)];
							//for(int q = 0; q < 8; q++)
							//	std::cout << verts[q]->Coords()[0] << " " << verts[q]->Coords()[1] << " " << verts[q]->Coords()[2] << " " << verts[q]->LocalID() << std::endl;
							Cell c = CreateCell(ElementArray<Node>(this,verts,verts+8),nvf,numnodes,6).first;
							if( !poro.empty() ) c->RealDF(tagporo) = poro[(i*dims[1]+j)*dims[2]+k];
							if( !perm.empty() )
							{
								Storage::real_array arr_perm = c->RealArrayDF(tagperm);
								arr_perm[0] = perm[3*((i*dims[1]+j)*dims[2]+k)+0];
								arr_perm[1] = perm[3*((i*dims[1]+j)*dims[2]+k)+1];
								arr_perm[2] = perm[3*((i*dims[1]+j)*dims[2]+k)+2];
							}
						}
			}
			else if( gtype == ECL_GTYPE_ZCORN )
			{
				//std::vector<HandleType> block_nodes(dims[0]*dims[1]*dims[2]*8);
				//Storage::integer block_zcorn[8];
			}
		}
		else
		if(LFile.find(".pvtk") != std::string::npos) //this is legacy parallel vtk
		{
			int state = 0, np, nf = 0;
			std::string::size_type attrl, ql, qr, l;
			//~ if( m_state == Mesh::Serial ) SetCommunicator(INMOST_MPI_COMM_WORLD);
			std::string str, tag, attrval, path = "";
			std::vector<std::string> files;
			std::fstream stream(File.c_str(),std::ios::in);
			l = File.find_last_of("/\\");
			if( l != std::string::npos )
				path = File.substr(0,l+1);
			while(stream.good())
			{
				getline(stream,str,'>');
				l = str.find_first_not_of(" \t\r\n");
				if( l != std::string::npos ) str.substr(l,str.size()-l).swap(str);
				else str.clear();
				if( str.empty() ) continue;
				l = str.find_first_of(' ');
				tag = str.substr(1,l-1);
				//std::cout << " tag name " << tag << std::endl;
				if( state == 0 )
				{
					if( tag == "File" )
					{
						attrl = str.find("version=\"pvtk-1.0\"",l);
						if( attrl == std::string::npos ) 
						{
							std::cout << __FILE__ << ":" << __LINE__ << " version information not found in " << File << std::endl;
							throw BadFile;
						}
						attrl = str.find("dataType=\"vtkUnstructuredGrid\"",l);
						if( attrl == std::string::npos ) 
						{
							std::cout << __FILE__ << ":" << __LINE__ << " data type information not found in " << File << std::endl;
							throw BadFile;
						}
						attrl = str.find("numberOfPieces");
						if( attrl == std::string::npos ) 
						{
							std::cout << __FILE__ << ":" << __LINE__ << " number of files not found in " << File << std::endl;
							throw BadFile;
						}
						ql = str.find_first_of('"',attrl);
						qr = str.find_first_of('"',ql+1);
						attrval = str.substr(ql+1,qr-ql-1);
						//std::cout << "number of peaces: " << attrval << " l: " << ql << " r: " << qr << std::endl;
						np = atoi(attrval.c_str());
						//std::cout << "np: " << np << std::endl;;
						files.resize(np);
						state = 1;
					}
					else 
					{
						std::cout << __FILE__ << ":" << __LINE__ << " unexpected xml tag " << tag << " in " << File << " expected File instead " << std::endl;
						throw BadFile;
					}
				}
				else if( state == 1 )
				{
					if( tag == "Piece")
					{
						attrl = str.find("fileName");
						ql = str.find_first_of('"',attrl);
						qr = str.find_first_of('"',ql+1);
						attrval = str.substr(ql+1,qr-ql-1);
						//std::cout << "file" << nf << ":" << attrval << std::endl;
						files[nf++] = attrval;
						if( nf == np ) state = 2;
					}
					else 
					{
						std::cout << __FILE__ << ":" << __LINE__ << " unexpected xml tag " << tag << " in " << File << " expected Piece instead " << std::endl;
						throw BadFile;
					}
				}
				else if( state == 2 )
				{
					if( tag != "/File" ) 
					{
						std::cout << __FILE__ << ":" << __LINE__ << " unexpected xml tag " << tag << " in " << File << " expected /File instead " << std::endl;
						throw BadFile;
					}
					else break;
				}
			}
			for(int i = GetProcessorRank(); i < static_cast<int>(files.size()); i+= GetProcessorsNumber()) 
			{
				//std::cout << "load " << i << ": " << path + files[i] << " by " << GetProcessorRank() << std::endl;
				Load(path + files[i]);
			}
			ResolveShared();
		}
		else if(LFile.find(".vtk") != std::string::npos) //this is legacy vtk
		{
			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;
			while( state != R_QUIT )
			{
				
				if( fgets(readline,2048,f) == NULL )
				{
					state = R_QUIT;
					continue;
				}
				if( readline[strlen(readline)-1] == '\n' )
					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:
					{
						//printf("Attached info: %s\n",readline);
						state = R_DATATYPE;
						break;
					}
					case R_DATATYPE:
					{
						if( !strncmp(readline,"BINARY",6) )
							binary = 1;
						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);
							if( verbosity > 0 ) printf("Reading attribute %s.\n",attrname);
							if( filled < 2 ) 
							{
								printf("%s:%d found %d arguments to SCALARS field, must be >= 2\nline:\n%s\n",__FILE__,__LINE__,filled,readline); 
								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 << 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 )
							{
								Tag attr = CreateTag(attrname,t,read_into_cell,read_into_cell & ESET,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newcells.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
								datatags.push_back(attr);
							}
							if( read_into == 1 )
							{
								Tag attr = CreateTag(attrname,t,NODE,NONE,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size())/250,1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newnodes.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
							}
							break;
						}
						else if( !strcmp(dataname,"COLOR_SCALARS") )
						{	
							char attrname[1024];
							filled = sscanf(readline," COLOR_SCALARS %s %d",attrname,&nentries);
							if( verbosity > 0 ) printf("Reading attribute %s.\n",attrname);
							if( read_into == 2 )
							{
								Tag attr = CreateTag(attrname,DATA_REAL,read_into_cell,read_into_cell & ESET,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
								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%100 == 0)
									{
										printf("data %3.1f%%\r",(it*report_pace)/(1.0*newcells.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
								datatags.push_back(attr);
							}
							if( read_into == 1 )
							{
								Tag attr = CreateTag(attrname,DATA_REAL,NODE,NONE,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size()/250),1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newnodes.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
							}
							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);
							if( verbosity > 0 ) printf("Reading attribute %s.\n",attrname);
							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 )
							{
								Tag attr = CreateTag(attrname,t,read_into_cell,read_into_cell & ESET,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newcells.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
								datatags.push_back(attr);
							}
							if( read_into == 1 )
							{
								Tag attr = CreateTag(attrname,t,NODE,NONE,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size()/250),1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newnodes.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
							}
							break;
						}
						else if( !strcmp(dataname,"TEXTURE_COORDINATES") )
						{
							char attrname[1024];
							char attrtype[1024];
							int nentries;
							DataType t = DATA_BULK;
							filled = sscanf(readline," TEXTURE_COORDINATES %s %d %s",attrname,&nentries,attrtype);
							if( filled != 3 ) 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 )
							{
								Tag attr = CreateTag(attrname,t,read_into_cell,read_into_cell & ESET,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newcells.size()/250),1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newcells.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
								datatags.push_back(attr);
							}
							if( read_into == 1 )
							{
								Tag attr = CreateTag(attrname,t,NODE,NONE,nentries);
								unsigned report_pace = std::max<unsigned>(static_cast<unsigned>(newnodes.size()/250),1);
								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)
									{
										printf("data %3.1f%%\r",(it*100.0)/(1.0*newnodes.size()));
										fflush(stdout);
									}
								}
								filled = fscanf(f,"\n");
							}
							break;
						}
						else if( !strcmp(dataname,"FIELD") )
						{
							int nfields = 0;
							filled = sscanf(readline," FIELD %*s %d",&nfields);
							if( filled != 1 ) throw BadFile;
							for(int i = 0; i < nfields; i++)
							{
								char attrname[1024];
								char attrtype[1024];
								unsigned int nentries, ntuples;
								DataType t = DATA_BULK;
								filled = fscanf(f," %s %u %u %s\n",attrname,&nentries,&ntuples,attrtype);
								if( filled != 4 ) 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( read_into == 2 )
								{
									Tag attr  = CreateTag(attrname,t,read_into_cell,read_into_cell & ESET,nentries);
									if( ntuples != newcells.size() ) printf("number of tuples in field is not equal to number of cells\n");
									unsigned report_pace = std::max<unsigned>(ntuples/250,1);
									for(unsigned int it = 0; it < ntuples; it++)
									{
										if( t == DATA_INTEGER )
										{
											if( newcells[it] != InvalidHandle() )
											{
												Storage::integer_array attrdata = IntegerArray(newcells[it],attr);
												for(unsigned int jt = 0; jt < nentries; jt++) {int temp; filled = fscanf(f," %d",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
											} else for(unsigned 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(unsigned int jt = 0; jt < nentries; jt++) {double temp; filled = fscanf(f," %lf",&temp); if(filled != 1 ) throw BadFile; attrdata[jt] = temp;}
											} else for(unsigned 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)
										{
											printf("data %3.1f%%\r",(it*100.0)/(1.0*ntuples));
											fflush(stdout);
										}
									}
									filled = fscanf(f,"\n");
									datatags.push_back(attr);
								}
								if( read_into == 1 )
								{
									Tag attr = CreateTag(attrname,t,NODE,NONE,nentries);
									if( ntuples != newnodes.size() ) printf("number of tuples in field is not equal to number of nodes\n");
									unsigned report_pace = std::max<unsigned>(ntuples/250,1);
									for(unsigned int it = 0; it < ntuples; it++)
									{
										if( t == DATA_INTEGER )
										{
											Storage::integer_array attrdata = IntegerArray(newnodes[it],attr);
											for(unsigned 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);