/*! \file IoModel.cpp
 * \brief  file containing the methods that will help in processing the input data coming 
 * into ISSM, from Matlab, or through a binary file opened for reading.
 */

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>

#include "./classes.h"
#include "../shared/io/io.h"
#include "../shared/shared.h"

/*FUNCTION IoModel::IoModel(){{{*/
IoModel::IoModel(){
	this->fid=NULL;
	this->data=NULL;
	this->independents=NULL;
	this->independent_objects=NULL;
	this->constants=NULL;

	this->my_elements=NULL;
	this->my_vertices=NULL;

	this->dim=-1;
	this->numberofvertices=-1;
	this->numberofelements=-1;
	this->numberoffaces=-1;
	this->numberofedges=-1;
	this->elements=NULL;
	this->faces=NULL;
	this->edges=NULL;
	this->elementtoedgeconnectivity      =NULL;
	this->singlenodetoelementconnectivity=NULL;
	this->numbernodetoelementconnectivity=NULL;

	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
/*FUNCTION IoModel::IoModel(FILE*  iomodel_handle){{{*/
IoModel::IoModel(FILE* iomodel_handle){

	/*First, keep track of the file handle: */
	this->fid=iomodel_handle;

	/*Check that Enums are Synchronized*/
	this->CheckEnumSync();

	/*Initialize data: */
	this->data=xNew<IssmDouble*>(MaximumNumberOfDefinitionsEnum);
	for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++) this->data[i]=NULL;

	/*If we are running in AD mode, we need to start the trace and declare our independent variables now, 
	 *and prevent them from being erased during successive calls to iomodel->FetchConstants, iomodel->FetchData and 
	 iomodel->DeleteData:*/
	this->DeclareIndependents();

	/*Initialize and read constants:*/
	this->constants=new Parameters();
	this->FetchConstants(); /*this routine goes through the input file, and fetches bools, ints, IssmDoubles and strings only, nothing memory intensive*/

	/*Initialize permanent data: */
	this->my_elements = NULL;
	this->my_vertices = NULL;

	FetchData(&this->dim,MeshDimensionEnum);
	FetchData(&this->numberofvertices,MeshNumberofverticesEnum);
	FetchData(&this->numberofelements,MeshNumberofelementsEnum);
	FetchData(&this->elements,NULL,NULL,MeshElementsEnum);
	this->faces                           = NULL;
	this->edges                           = NULL;
	this->elementtoedgeconnectivity       = NULL;
	this->singlenodetoelementconnectivity = NULL;
	this->numbernodetoelementconnectivity = NULL;

	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
/*FUNCTION IoModel::~IoModel(){{{*/
IoModel::~IoModel(){

	/*Some checks in debugging mode*/
	#ifdef _ISSM_DEBUG_
	if(this->data){
		for(int i=0;i<MaximumNumberOfDefinitionsEnum;i++){
			if(this->data[i]){
				_printf0_("Info: previous pointer of " << EnumToStringx(i) << " has not been freed (DeleteData has not been called)\n");
			}
		}
	}
	#endif

	if(this->constants) delete this->constants;

	xDelete<IssmDouble*>(this->data);
	xDelete<bool>(this->independents);
	if(this->independent_objects)delete this->independent_objects;

	xDelete<bool>(this->my_elements);
	xDelete<int>(this->my_vertices);

	xDelete<int>(this->elements);
	xDelete<int>(this->faces);
	xDelete<int>(this->edges);
	xDelete<int>(this->elementtoedgeconnectivity);
	xDelete<int>(this->singlenodetoelementconnectivity);
	xDelete<int>(this->numbernodetoelementconnectivity);
}
/*}}}*/

/*FUNCTION IoModel::CheckEnumSync{{{*/
void  IoModel::CheckEnumSync(void){

	int my_rank;
	int record_enum = 0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Check that some fields have been allocated*/
	_assert_(this->fid || my_rank);

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){ //cpu 0

		/*First set FILE* position to the beginning of the file: */
		fseek(this->fid,0,SEEK_SET);

		/*Get first Enum*/
		if(fread(&record_enum,sizeof(int),1,this->fid)==0){
			_error_("Marshalled file is empty");
		}
		else{
			if(record_enum!=MaximumNumberOfDefinitionsEnum){
				_printf0_("\n");
				_printf0_("=========================================================================\n");
				_printf0_(" Enums in marshalled file are not compatible with compiled code          \n");
				_printf0_("                                                                         \n");
				_printf0_("   * If you are running ISSM on a remote cluster:                        \n");
				_printf0_("     make sure that you are using the same version of ISSM on your local \n");
				_printf0_("     machine and remote cluster (you might need to run svn update)       \n");
				_printf0_("   * If you are running ISSM on your local machine:                      \n");
				_printf0_("     make sure that all the code is compiled (modules and executables)   \n");
				_printf0_("   * If you are a developer and just added a new Enum:                   \n");
				_printf0_("     you might need to run ./Synchronize.sh in src/c/EnumDefinitions     \n");
				_printf0_("     and recompile                                                       \n");
				_printf0_("=========================================================================\n");
				_printf0_("\n");
				_error_("Enums not consistent (See error message above)");
			}
		}
	}
}
/*}}}*/
/*FUNCTION IoModel::Constant(bool* poutput,int constant_enum){{{*/
void IoModel::Constant(bool* poutput,int constant_enum){

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
/*FUNCTION IoModel::Constant(int* poutput,int constant_enum){{{*/
void IoModel::Constant(int* poutput,int constant_enum){

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
/*FUNCTION IoModel::Constant(IssmDouble* poutput,int constant_enum){{{*/
void IoModel::Constant(IssmDouble* poutput,int constant_enum){

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
/*FUNCTION IoModel::Constant(char** poutput,int constant_enum){{{*/
void IoModel::Constant(char** poutput,int constant_enum){

	_assert_(constant_enum>=0);
	_assert_(this->constants);

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
/*FUNCTION IoModel::CopyConstantObject{{{*/
Param* IoModel::CopyConstantObject(int constant_enum){

	_assert_(this->constants);

	/*Find constant*/
	Param* param=dynamic_cast<Param*>(this->constants->FindParamObject(constant_enum));
	if(!param) _error_("Constant " << EnumToStringx(constant_enum) << " not found in iomodel");

	return dynamic_cast<Param*>(param->copy());
}
/*}}}*/
/*FUNCTION IoModel::Data{{{*/
IssmDouble* IoModel::Data(int data_enum){

	_assert_(data_enum<MaximumNumberOfDefinitionsEnum);
	_assert_(data_enum>=0);

	return this->data[data_enum];
}
/*}}}*/
/*FUNCTION IoModel::DeclareIndependents{{{*/
void IoModel::DeclareIndependents(void){

	int  i;
	bool autodiff = false;
	int  num_independent_objects;

	int *names = NULL;
	int *types = NULL;

	int  dummy;
	bool keep=false;

	int my_rank=IssmComm::GetRank();

	/*Initialize array detecting whether data[i] is an independent AD mode variable: */
	this->independents=xNew<bool>(MaximumNumberOfDefinitionsEnum);
	for(i=0;i<MaximumNumberOfDefinitionsEnum;i++) this->independents[i]=false;

	this->FetchData(&autodiff,AutodiffIsautodiffEnum);
	if(autodiff){

		#ifdef _HAVE_ADOLC_
		/*Start trace*/
		this->FetchData(&keep,AutodiffKeepEnum);
		if(keep)trace_on(my_rank,1);
		else    trace_on(my_rank);

		/*build dataset made of independent objects:*/
		this->independent_objects=new DataSet();
		this->FetchData(&num_independent_objects,AutodiffNumIndependentObjectsEnum);
		if(num_independent_objects){
			this->FetchData(&names,&dummy,&dummy,AutodiffIndependentObjectNamesEnum);
			this->FetchData(&types,&dummy,&dummy,AutodiffIndependentObjectTypesEnum);

			/*create independent objects, and at the same time, fetch the corresponding independent variables, 
			 *and declare them as such in ADOLC: */
			for(i=0;i<num_independent_objects;i++){

				IndependentObject* independent_object=NULL;
				independent_object=new IndependentObject(names[i],types[i]);

				/*add to independent_objects dataset:*/
				this->independent_objects->AddObject(independent_object);

				/*now go fetch the independent variable: */
				independent_object->FetchIndependent(this); //supply the pointer to iomodel.
			}
			xDelete<int>(names);
			xDelete<int>(types);
		}
		#else
		/*if we asked for AD computations, we have a problem!: */
		_error_("Cannot carry out AD mode computations without support of ADOLC compiled in!");
		#endif
	}
	else this->independent_objects=NULL;

}
/*}}}*/
/*FUNCTION IoModel::DeleteData(int num,...){{{*/
void  IoModel::DeleteData(int num,...){

	va_list ap;
	int     dataenum;
	int     i;
	DoubleMatParam* parameter=NULL;

	/*Go through the entire list of enums and delete the corresponding data from the iomodel-data dataset: */

	va_start(ap,num);
	for(i = 0; i <num; i++){
		dataenum=va_arg(ap, int);
		_assert_(dataenum<MaximumNumberOfDefinitionsEnum);

		/*do not erase independent variables for the AD mode computations!: */
		if (!this->independents[dataenum]) xDelete<IssmDouble>(this->data[dataenum]);
	}
	va_end(ap);
} /*}}}*/
/*FUNCTION IoModel::DeleteData(IssmDouble* {{{*/
void  IoModel::DeleteData(IssmDouble* vector, int dataenum){

	/*do not erase independent variables for the AD mode computations!: */
	if(vector)if (!this->independents[dataenum]) xDelete<IssmDouble>(vector);

} /*}}}*/
/*FUNCTION IoModel::FetchConstants{{{*/
void  IoModel::FetchConstants(void){

	int my_rank;

	/*record descriptions; */
	int record_enum;
	int record_length;
	int record_code; //1 to 7 number

	/*records: */
	int  booleanint=0;
	int  integer=0;
	IssmPDouble pscalar=0;
	IssmDouble scalar=0;
	char* string=NULL;
	int   string_size;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Check that some fields have been allocated*/
	_assert_(this->fid || my_rank);
	_assert_(this->constants);

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){ //cpu 0{{{

		/*First set FILE* position to the beginning of the file: */
		fseek(this->fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifiers (bool,int,IssmDouble or string): */
		for(;;){
			if(fread(&record_enum,sizeof(int),1,this->fid)==0){

				/*Ok, we have reached the end of the file. break: */
				record_code=0; //0 means bailout
				ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm());  /*tell others cpus we are bailing: */
				break;
			}
			else{

				/* Read the record length and the data type code: */
				if(fread(&record_length,sizeof(int),1,this->fid)!=1) _error_("Cound not read record_length");
				if(fread(&record_code  ,sizeof(int),1,this->fid)!=1) _error_("Cound not read record_code");

				/*Tell other cpus what we are doing: */
				ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm());  /*tell other cpus what we are going to do: */

				/*Tell other cpus the name of the data, then branch according to the data type: */
				ISSM_MPI_Bcast(&record_enum,1,ISSM_MPI_INT,0,IssmComm::GetComm());  
				ISSM_MPI_Bcast(&record_length,1,ISSM_MPI_INT,0,IssmComm::GetComm());  

				switch(record_code){
					case 1: 
						/*Read the boolean and broadcast it to other cpus:*/
						if(fread(&booleanint,sizeof(int),1,this->fid)!=1) _error_("could not read boolean ");
						ISSM_MPI_Bcast(&booleanint,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

						/*create BoolParam: */
						this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to boolean

						break;
					case 2:
						/*Read the integer and broadcast it to other cpus:*/
						if(fread(&integer,sizeof(int),1,this->fid)!=1) _error_("could not read integer ");
						ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

						/*create IntParam: */
						this->constants->AddObject(new IntParam(record_enum,integer));

						break;
					case 3:
						/*Read the scalar and broadcast it to other cpus. However, if this record has already 
						 * been read in DeclareIndependents, we should not re-read it, but grab it from the iomodel->data 
						 * slots: */
						if(this->independents[record_enum]){
							scalar=*(this->data[record_enum]);
						}
						else{
							if(fread(&pscalar,sizeof(IssmPDouble),1,this->fid)!=1) _error_("could not read scalar ");
							ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
							scalar=pscalar;
						}

						/*create DoubleParam: */
						this->constants->AddObject(new DoubleParam(record_enum,scalar));

						break;
					case 4: 
						/*We have to read a string from disk. First read the dimensions of the string, then the string: */
						if(fread(&string_size,sizeof(int),1,this->fid)!=1) _error_("could not read length of string ");
						ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

						if(string_size){
							string=xNew<char>(string_size+1);
							string[string_size]='\0';

							/*Read string, then broadcast: */
							if(fread(string,string_size*sizeof(char),1,this->fid)!=1)_error_(" could not read string ");
							ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
						}
						else{
							string=xNew<char>(1);
							string[0]='\0';
						}

						/*Add string to parameters: */
						this->constants->AddObject(new StringParam(record_enum,string));

						/*Free string*/
						xDelete<char>(string);

						break;
					case 5: 
							/*We are not interested in this record, too memory intensive. Skip it: */
							/*skip: */
							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
							fseek(fid,record_length,SEEK_CUR);
							break;
					case 6: 
							/*We are not interested in this record, too memory intensive. Skip it: */
							/*skip: */
							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
							fseek(fid,record_length,SEEK_CUR);
							break;
					case 7: 
							/*We are not interested in this record, too memory intensive. Skip it: */
							/*skip: */
							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
							fseek(fid,record_length,SEEK_CUR);
							break;

					case 8: 
							/*We are not interested in this record, too memory intensive. Skip it: */
							/*skip: */
							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
							fseek(fid,record_length,SEEK_CUR);
							break;

					case 9: 
							/*We are not interested in this record, too memory intensive. Skip it: */
							/*skip: */
							fseek(fid,-sizeof(int),SEEK_CUR); //backtrak 1 integer
							fseek(fid,record_length,SEEK_CUR);
							break;

					default: 
						_error_("unknown record type:" << record_code); 
						break;;
				}
			}
		}
	} //}}}
	else{ //cpu ~0 {{{
		for(;;){ //wait on cpu 0
			ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm());  /*get from cpu 0 what we are going to do: */
			if(record_code==0){
				break; //we are done, break from the loop
			}
			else{
				ISSM_MPI_Bcast(&record_enum,1,ISSM_MPI_INT,0,IssmComm::GetComm());   //get from cpu 0 name of the data
				ISSM_MPI_Bcast(&record_length,1,ISSM_MPI_INT,0,IssmComm::GetComm());  
				switch(record_code){
				case 1: 
					/*boolean. get it from cpu 0 */
					ISSM_MPI_Bcast(&booleanint,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

					/*create BoolParam: */
					this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean
					break;

				case 2:
					/*integer. get it from cpu 0 */
					ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

					/*create IntParam: */
					this->constants->AddObject(new IntParam(record_enum,integer));

					break;
				case 3:
					/*scalar. get it from cpu 0 */
					ISSM_MPI_Bcast(&pscalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
					scalar=pscalar;
					/*create DoubleParam: */
					this->constants->AddObject(new DoubleParam(record_enum,scalar));

					break;
				case 4: 
					ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
					if(string_size){
						string=xNew<char>((string_size+1));
						string[string_size]='\0';

						/*Read string from cpu 0: */
						ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
					}
					else{
						string=xNew<char>(1);
						string[0]='\0';
					}
					/*Add string to parameters: */
					this->constants->AddObject(new StringParam(record_enum,string));

					/*Free string*/
					xDelete<char>(string);

					break;
				case 5: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 6: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 7: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 8: break; //do nothing. not interested in this type of data, which is memory intensive.
				case 9: break; //do nothing. not interested in this type of data, which is memory intensive.

				default: 
					_error_("unknown record type:" << record_code); 
					break;;
				}

			}
		}
	} //}}}
}
/*}}}*/
/*FUNCTION IoModel::FetchData(bool*     pbool,int data_enum){{{*/
void  IoModel::FetchData(bool* pboolean,int data_enum){

	int my_rank;

	/*output: */
	int   booleanint;
	int   code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=1)_error_("expecting a boolean for enum " << EnumToStringx(data_enum));

	/*We have to read a boolean from disk. */
	if(my_rank==0){  
		if(fread(&booleanint,sizeof(int),1,fid)!=1) _error_("could not read boolean ");
	}
	ISSM_MPI_Bcast(&booleanint,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*cast to bool: */
	/*Assign output pointers: */
	*pboolean=(bool)booleanint;

}
/*}}}*/
/*FUNCTION IoModel::FetchData(int*      pinteger,int data_enum){{{*/
void  IoModel::FetchData(int* pinteger,int data_enum){

	int my_rank;

	/*output: */
	int   integer;
	int   code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=2)_error_("expecting an integer for enum " << EnumToStringx(data_enum));

	/*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */
	if(my_rank==0){  
		if(fread(&integer,sizeof(int),1,fid)!=1) _error_("could not read integer ");
	}

	ISSM_MPI_Bcast(&integer,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers: */
	*pinteger=integer;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(IssmDouble*   pscalar,int data_enum){{{*/
void  IoModel::FetchData(IssmDouble* pscalar,int data_enum){

	int my_rank;

	/*output: */
	IssmPDouble   scalar;
	int      code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=3)_error_("expecting a IssmDouble for enum " << EnumToStringx(data_enum));

	/*We have to read a scalar from disk. First read the dimensions of the scalar, then the scalar: */
	if(my_rank==0){
		if(fread(&scalar,sizeof(IssmPDouble),1,fid)!=1)_error_("could not read scalar ");
	}
	ISSM_MPI_Bcast(&scalar,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 

	/*Assign output pointers: */
	*pscalar=scalar;

}
/*}}}*/
/*FUNCTION IoModel::FetchData(char**    pstring,int data_enum){{{*/
void  IoModel::FetchData(char** pstring,int data_enum){

	int my_rank;

	/*output: */
	char* string=NULL;
	int   string_size;
	int code=0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=4)_error_("expecting a string for enum " << EnumToStringx(data_enum));

	/*Now fetch: */

	/*We have to read a string from disk. First read the dimensions of the string, then the string: */
	if(my_rank==0){  
		if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
	}

	ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Now allocate string: */
	if(string_size){
		string=xNew<char>((string_size+1));
		string[string_size]='\0';

		/*Read string on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(string,string_size*sizeof(char),1,fid)!=1)_error_(" could not read string ");
		}
		ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
	}
	else{
		string=xNew<char>(1);
		string[0]='\0';
	}

	/*Assign output pointers: */
	*pstring=string;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(int**     pintegerematrix,int* pM,int* pN,int data_enum){{{*/
void  IoModel::FetchData(int** pmatrix,int* pM,int* pN,int data_enum){

	int my_rank;
	int i,j;

	/*output: */
	int M,N;
	IssmPDouble* matrix=NULL;
	int*    integer_matrix=NULL;
	int code=0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum));

	/*Now fetch: */

	/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
	/*numberofelements: */
	if(my_rank==0){  
		if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
	}

	ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	if(my_rank==0){  
		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
	}
	ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());

	/*Now allocate matrix: */
	if(M*N){
		matrix=xNew<IssmPDouble>(M*N);

		/*Read matrix on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
		}

		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
	}

	/*Now cast to integer: */
	if(M*N){
		integer_matrix=xNew<int>(M*N);
		for (i=0;i<M;i++){
			for (j=0;j<N;j++){
				integer_matrix[i*N+j]=(int)matrix[i*N+j];
			}
		}
	}
	else{
		integer_matrix=NULL;
	}
	/*Free ressources:*/
	xDelete<IssmPDouble>(matrix);

	/*Assign output pointers: */
	*pmatrix=integer_matrix;
	if (pM)*pM=M;
	if (pN)*pN=N;

}
/*}}}*/
/*FUNCTION IoModel::FetchData(IssmDouble**  pIssmDoublematrix,int* pM,int* pN,int data_enum){{{*/
void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN,int data_enum){

	int my_rank;

	/*output: */
	int M,N;
	IssmPDouble* matrix=NULL;
	int code=0;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	if((code!=5) && (code!=6) && (code!=7))_error_("expecting a IssmDouble, integer or boolean matrix for enum " << EnumToStringx(data_enum));

	/*Now fetch: */

	/*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
	/*numberofelements: */
	if(my_rank==0){  
		if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
	}
	ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	if(my_rank==0){  
		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
	}
	ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Now allocate matrix: */
	if(M*N){
		matrix=xNew<IssmPDouble>(M*N);

		/*Read matrix on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
		}
		ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 

		_assert_(this->independents);
		if (this->independents[data_enum]){
			/*this data has already been checked out! So cancel all that we've done here, and return 
			 * the data[data_enum] directly: */
			*pmatrix=this->data[data_enum];
		}
		else{
			*pmatrix=xNew<IssmDouble>(M*N);
			for (int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
		}
		xDelete<IssmPDouble>(matrix);
	}
	else
	  *pmatrix=NULL;
	/*Assign output pointers: */
	if (pM)*pM=M;
	if (pN)*pN=N;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(char***   pstrings,int* pnumstrings,int data_enum){{{*/
void  IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){

	int my_rank;

	int i;

	/*output: */
	int   numstrings=0;
	char** strings=NULL;

	/*intermediary: */
	char* string=NULL;
	int   string_size;
	int   code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);

	if(code!=9)_error_("expecting a string array for enum " << EnumToStringx(data_enum));

	/*We have to read a bunch of strings from disk. First read the number of strings, and allocate: */
	if(my_rank==0){  
		if(fread(&numstrings,sizeof(int),1,fid)!=1) _error_("could not read length of string array");
	}
	ISSM_MPI_Bcast(&numstrings,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Now allocate string array: */
	if(numstrings){
		strings=xNew<char*>(numstrings);
		for(i=0;i<numstrings;i++)strings[i]=NULL;

		/*Go through strings, and read: */
		for(i=0;i<numstrings;i++){

			if(my_rank==0){  
				if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
			}
			ISSM_MPI_Bcast(&string_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
			if(string_size){
				string=xNew<char>((string_size+1));
				string[string_size]='\0';

				/*Read string on node 0, then broadcast: */
				if(my_rank==0){  
					if(fread(string,string_size*sizeof(char),1,fid)!=1)_error_(" could not read string ");
				}
				ISSM_MPI_Bcast(string,string_size,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
			}
			else{
				string=xNew<char>(1);
				string[0]='\0';
			}

			strings[i]=string;
		}
	}

	/*Assign output pointers: */
	*pstrings=strings;
	*pnumstrings=numstrings;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pM,int data_enum){{{*/
void  IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,int data_enum){

	int i;

	int my_rank;

	/*output: */
	IssmDouble** matrices=NULL;
	int*     mdims=NULL;
	int*     ndims=NULL;
	int      numrecords=0;

	/*intermediary: */
	int     M, N;
	IssmPDouble *matrix = NULL;
	int     code;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	if(code!=8)_error_("expecting a IssmDouble mat array for enum " << EnumToStringx(data_enum));

	/*Now fetch: */
	if(my_rank==0){  
		if(fread(&numrecords,sizeof(int),1,fid)!=1) _error_("could not read number of records in matrix array ");
	}
	ISSM_MPI_Bcast(&numrecords,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	if(numrecords){

		/*Allocate matrices :*/
		matrices=xNew<IssmDouble*>(numrecords);
		mdims=xNew<int>(numrecords);
		ndims=xNew<int>(numrecords);

		for(i=0;i<numrecords;i++){
			matrices[i]=NULL;
			mdims[i]=0;
			ndims[i]=0;
		}

		/*Loop through records and fetch matrix: */
		for(i=0;i<numrecords;i++){

			if(my_rank==0){  
				if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
			}
			ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

			if(my_rank==0){  
				if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns in " << i << "th matrix of matrix array");
			}
			ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

			/*Now allocate matrix: */
			if(M*N){
				matrix=xNew<IssmPDouble>(M*N);

				/*Read matrix on node 0, then broadcast: */
				if(my_rank==0){  
					if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
				}

				ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
				matrices[i]=xNew<IssmDouble>(M*N);
				for (int j=0;j<M*N;++j) {matrices[i][j]=matrix[j];}
				xDelete<IssmPDouble>(matrix);
			}
			else
			  matrices[i]=NULL;
			/*Assign: */
			mdims[i]=M;
			ndims[i]=N;
		}
	}

	/*Assign output pointers: */
	*pmatrices=matrices;
	*pmdims=mdims;
	*pndims=ndims;
	*pnumrecords=numrecords;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(Option**  poption,int data_enum){{{*/
void  IoModel::FetchData(Option** poption,int index){

	/*output: */
	int     code;
	char   *name        = NULL;

	/*First get option name*/
	this->FetchData(&name,index);

	/*Get option value*/
	fid=this->SetFilePointerToData(&code,NULL,index+1);
	switch(code){
		case 3: {//IssmDouble
					GenericOption<IssmDouble>* option;
					IssmDouble value;
					FetchData(&value,index+1);
					option = new GenericOption<IssmDouble>();
					option->value = value;
					option->name  = name;
					option->numel = 1;
					option->ndims = 1;
					option->size  = NULL;
					/*Assign output pointers: */
					*poption=option;
					break;
				}
		case 4: {//char
					GenericOption<char*>* option;
					char* value = NULL;
					FetchData(&value,index+1);
					option = new GenericOption<char*>();
					option->value = value;
					option->name  = name;
					option->numel = 1;
					option->ndims = 1;
					option->size  = NULL;
					*poption=option;
					break;
				}
		default:
			  _error_("Option of format " << code << " not supported yet");
	}

}
/*}}}*/
/*FUNCTION IoModel::FetchData(int num,...){{{*/
void  IoModel::FetchData(int num,...){

	va_list     ap;
	int         dataenum;
	IssmDouble *matrix   = NULL;
	int         M,N;
	int         i;

	/*Go through the entire list of enums and fetch the corresponding data. Add it to the iomodel->data dataset. Everything
	 *we fetch is a IssmDouble* : */

	va_start(ap,num);
	for(i=0; i<num; i++){

		dataenum=va_arg(ap, int);
		_assert_(dataenum<MaximumNumberOfDefinitionsEnum); 
		_assert_(dataenum>=0);

		if (this->independents[dataenum]){
			/*this data has already been checked out! Continue: */
			continue;
		}

		/*Some checks in debugging mode*/
		/*{{{*/
		#ifdef _ISSM_DEBUG_
		_assert_(dataenum<MaximumNumberOfDefinitionsEnum);
		if(this->data[dataenum]){
			_error_("Info: trying to fetch " << EnumToStringx(dataenum) << " but previous pointer has not been freed (DeleteData has not been called)");
		}
		#endif
		/*}}}*/

		/*Add to this->data: */
		this->FetchData(&matrix,&M,&N,dataenum);
		this->data[dataenum]=matrix;
	}
	va_end(ap);

}
/*}}}*/
/*FUNCTION IoModel::FetchDataToInput{{{*/
void IoModel::FetchDataToInput(Elements* elements,int vector_enum,int default_vector_enum,IssmDouble default_value){

	/*intermediary: */
	int     i;
	bool    defaulting    = false;
	bool    transient     = false;

	FILE   *fid           = NULL;
	int     code          = 0;
	int     vector_layout = 0;
	int     counter;
	int     numberofelements;

	/*variables being fetched: */
	bool        boolean;
	int         integer;
	IssmDouble  scalar;
	char       *string           = NULL;
	IssmDouble *IssmDoublevector = NULL;
	int         M,N;

	/*Fetch parameters: */
	this->constants->FindParam(&numberofelements,MeshNumberofelementsEnum);

	/*First of, find the record for the enum, and get code  of data type: */
	fid=this->SetFilePointerToData(&code, &vector_layout,vector_enum);

	switch(code){
		case 1: //boolean constant.  {{{
			this->FetchData(&boolean,vector_enum);

			/*Add boolean constant input to all elements: */
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(counter));
					element->InputCreate((IssmDouble)boolean,vector_enum,code);
					counter++;
				}
			}
			break; /*}}}*/
		case 2: //integer constant.  {{{
			this->FetchData(&integer,vector_enum);

			/*Add integer constant input to all elements: */
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(counter));
					element->InputCreate((IssmDouble)integer,vector_enum,code);
					counter++;
				}
			}
			break; /*}}}*/
		case 3: //IssmDouble constant.  {{{
			this->FetchData(&scalar,vector_enum);

			/*Add IssmDouble constant input to all elements: */
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(counter));
					element->InputCreate(scalar,vector_enum,code); 
					counter++;
				}
			}
			break; /*}}}*/
		case 5: //boolean vector.  {{{

			/*Fetch vector:*/
			this->FetchData(&IssmDoublevector,&M,&N,vector_enum); //we still have a IssmDoublevector, because it might include times in transient mode
			/*Check we got something, otherwise fetch default: */
			if(IssmDoublevector){
				defaulting=false;  //we are not defaulting, because  we do have a vector
			}
			else{
				/*do we have a default enum to go fetch? */
				if(default_vector_enum!=NoneEnum){
					/*yes. fetch it: */
					this->FetchData(&IssmDoublevector,&M,&N,default_vector_enum);
					if(IssmDoublevector){
						defaulting=false;  //we are not defaulting, because  we do have a vector
					}
					else{
						/*even the default vector is non existent. we are defaulting to the default value: */
						defaulting=true;
					}
				}
				else{
					/*we are left with the default value: */
					defaulting=true;
				}
			}

			/*Create inputs:*/
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(counter));
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(IssmDoublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}
			break; /*}}}*/
		case 6: //int vector{{{

			/*Fetch vector:*/
			this->FetchData(&IssmDoublevector,&M,&N,vector_enum); //we still have a IssmDoublevector, because it might include times in transient mode
			/*Check we got something, otherwise fetch default: */
			if(IssmDoublevector){
				defaulting=false;  //we are not defaulting, because  we do have a vector
			}
			else{
				/*do we have a default enum to go fetch? */
				if(default_vector_enum!=NoneEnum){
					/*yes. fetch it: */
					this->FetchData(&IssmDoublevector,&M,&N,default_vector_enum);
					if(IssmDoublevector){
						defaulting=false;  //we are not defaulting, because  we do have a vector
					}
					else{
						/*even the default vector is non existent. we are defaulting to the default value: */
						defaulting=true;
					}
				}
				else{
					/*we are left with the default value: */
					defaulting=true;
				}
			}

			/*Create inputs:*/
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(counter));
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(IssmDoublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}
			break; /*}}}*/
		case 7: //IssmDouble vector{{{

			/*Fetch vector:*/
			this->FetchData(&IssmDoublevector,&M,&N,vector_enum);
			/*Check we got something, otherwise fetch default: */
			if(IssmDoublevector){
				defaulting=false;  //we are not defaulting, because  we do have a vector
			}
			else{
				/*do we have a default enum to go fetch? */
				if(default_vector_enum!=NoneEnum){
					/*yes. fetch it: */
					this->FetchData(&IssmDoublevector,&M,&N,default_vector_enum);
					if(IssmDoublevector){
						defaulting=false;  //we are not defaulting, because  we do have a vector
					}
					else{
						/*even the default vector is non existent. we are defaulting to the default value: */
						defaulting=true;
					}
				}
				else{
					/*we are left with the default value: */
					defaulting=true;
				}
			}

			/*Create inputs:*/
			counter=0;
			for (i=0;i<numberofelements;i++){
				if(this->my_elements[i]){
					Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(counter));
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(IssmDoublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}

			break; /*}}}*/
		default: /*{{{*/
			_error_("data code " << code << " not supported yet!");
			break;
			/*}}}*/
	}
	/*Free ressources. Pay attention to not freeing an AD mode independent variable though!:*/
	if (!this->independents[vector_enum] && !this->independents[default_vector_enum]) xDelete<IssmDouble>(IssmDoublevector);
	xDelete<char>(string);
}
/*}}}*/
/*FUNCTION IoModel::LastIndex{{{*/
void IoModel::LastIndex(int *pindex){

	int my_rank;
	int lastindex,index;
	int record_length;

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){

		/*First set FILE* position to the beginning of the file: */
		fseek(fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifier: */
		for(;;){
			/*Read enum for this size of first string name: */
			if(fread(&index,sizeof(int),1,fid)==0){
				/*Ok, we have reached the end of the file. break: */
				break;
			}

			/*read the record length, and use it to skip this record: */
			if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
			fseek(fid,record_length,SEEK_CUR);
			lastindex=index;
		}
	}
	/*Broadcast code and vector type: */
	ISSM_MPI_Bcast(&lastindex,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers:*/
	*pindex=lastindex;
}
/*}}}*/
/*FUNCTION IoModel::SetFilePointerToData{{{*/
FILE* IoModel::SetFilePointerToData(int* pcode,int* pvector_type, int data_enum){

	int my_rank;

	int found         = 0;
	int record_enum;
	int record_length;
	int record_code;       //1 to 7 number
	int vector_type   = 0; //nodal or elementary

	/*recover my_rank:*/
	my_rank=IssmComm::GetRank();

	/*Go find in the binary file, the position of the data we want to fetch: */
	if(my_rank==0){

		/*First set FILE* position to the beginning of the file: */
		fseek(fid,0,SEEK_SET);

		/*Now march through file looking for the correct data identifier: */
		for(;;){
			/*Read enum for this size of first string name: */
			if(fread(&record_enum,sizeof(int),1,fid)==0){
				/*Ok, we have reached the end of the file. break: */
				found=0;
				break;
			}

			/*Is this the record sought for? : */
			if (data_enum==record_enum){
				/*Ok, we have found the correct string. Pass the record length, and read data type code: */
				fseek(fid,sizeof(int),SEEK_CUR);
				if(fread(&record_code,sizeof(int),1,fid)!=1) _error_("Could not read record_code");

				/*if record_code points to a vector, get its type (nodal or elementary): */
				if(5<=record_code && record_code<=7){
					if(fread(&vector_type,sizeof(int),1,fid)!=1) _error_("Could not read vector_type");
				}
				found=1;
				break;
			}
			else{
				/*This is not the correct string, read the record length, and use it to skip this record: */
				if(fread(&record_length,sizeof(int),1,fid)!=1) _error_("Could not read record_length");
				/*skip: */
				fseek(fid,record_length,SEEK_CUR);
			}
		}
	}
	ISSM_MPI_Bcast(&found,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
	if(!found)_error_("could not find data with name " << EnumToStringx(data_enum) << " in binary file");

	/*Broadcast code and vector type: */
	ISSM_MPI_Bcast(&record_code,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
	ISSM_MPI_Bcast(&vector_type,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
	if(record_code==5) ISSM_MPI_Bcast(&vector_type,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 

	/*Assign output pointers:*/
	*pcode=record_code;
	if(pvector_type)*pvector_type=vector_type;

	return fid;
}
/*}}}*/
