/*! \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 "./objects.h"
#include "../io/io.h"
#include "./Container/Parameters.h"
#include "../shared/shared.h"
#include "../io/io.h"
#include "../include/include.h"

/*FUNCTION IoModel::IoModel(){{{1*/
IoModel::IoModel(){
	this->fid=NULL;
	this->data=NULL;
	this->constants=NULL;
	
	this->my_elements=NULL;
	this->my_nodes=NULL;
	this->my_vertices=NULL;
	this->singlenodetoelementconnectivity=NULL;
	this->numbernodetoelementconnectivity=NULL;
	
	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
/*FUNCTION IoModel::IoModel(FILE*  iomodel_handle){{{1*/
IoModel::IoModel(FILE* iomodel_handle){
	
	/*First, keep track of the file handle: */
	this->fid=iomodel_handle;

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

	/*Initialize data: */
	this->data=(double**)xmalloc(MaximumNumberOfEnums*sizeof(double*));
	for(int i=0;i<MaximumNumberOfEnums;i++) this->data[i]=NULL;
	
	/*Initialize permanent data: */
	this->my_elements=NULL;
	this->my_nodes=NULL;
	this->my_vertices=NULL;
	this->singlenodetoelementconnectivity=NULL;
	this->numbernodetoelementconnectivity=NULL;
	
	this->nodecounter=0;
	this->loadcounter=0;
	this->constraintcounter=0;
}
/*}}}*/
/*FUNCTION IoModel::~IoModel(){{{1*/
IoModel::~IoModel(){

	delete this->constants;

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

	xfree((void**)&this->data);

	xfree((void**)&this->my_elements);
	xfree((void**)&this->my_nodes);
	xfree((void**)&this->my_vertices);
	xfree((void**)&this->singlenodetoelementconnectivity);
	xfree((void**)&this->numbernodetoelementconnectivity);
}
/*}}}*/

/*FUNCTION IoModel::Constant(bool* poutput,int constant_enum){{{1*/
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){{{1*/
void IoModel::Constant(int* poutput,int constant_enum){

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

	this->constants->FindParam(poutput,constant_enum);
}
/*}}}*/
/*FUNCTION IoModel::Constant(double* poutput,int constant_enum){{{1*/
void IoModel::Constant(double* 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){{{1*/
void IoModel::Constant(char** poutput,int constant_enum){

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

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

	_assert_(this->constants);

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

	return (Param*)param->copy();
}
/*}}}*/
/*FUNCTION IoModel::Data{{{1*/
double* IoModel::Data(int data_enum){

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

	return this->data[data_enum];
}
/*}}}*/
/*FUNCTION IoModel::DeleteData{{{1*/
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<MaximumNumberOfEnums);
		xfree((void**)&this->data[dataenum]);
	}
	va_end(ap);
} /*}}}*/
/*FUNCTION IoModel::FetchConstants{{{1*/
void  IoModel::FetchConstants(void){

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

	/*records: */
	int  booleanint=0;
	int  integer=0;
	double scalar=0;
	char* string=NULL;
	int   string_size;

	/*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{{{2
	
		/*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,double 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
				#ifdef _HAVE_MPI_
				MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD);  /*tell others cpus we are bailing: */
				#endif
				break;
			}
			else{
			
				/* Read the record length and the data type code: */
				fread(&record_length,sizeof(int),1,this->fid);
				fread(&record_code,sizeof(int),1,this->fid);
					
				#ifdef _HAVE_MPI_
				/*Tell other cpus what we are doing: */
				MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD);  /*tell other cpus what we are going to do: */

				/*Tell other cpus the name of the data, then branch according to the data type: */
				MPI_Bcast(&record_enum,1,MPI_INT,0,MPI_COMM_WORLD);  
				MPI_Bcast(&record_length,1,MPI_INT,0,MPI_COMM_WORLD);  
				#endif
				
				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 ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&booleanint,1,MPI_INT,0,MPI_COMM_WORLD); 
						#endif

						/*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 ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&integer,1,MPI_INT,0,MPI_COMM_WORLD); 
						#endif

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

						break;
					case 3:
						/*Read the scalar and broadcast it to other cpus:*/
						if(fread(&scalar,sizeof(double),1,this->fid)!=1) _error_(" could not read scalar ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
						#endif

						/*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 ");
						#ifdef _HAVE_MPI_
						MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
						#endif

						if(string_size){
							string=(char*)xmalloc((string_size+1)*sizeof(char));
							string[string_size]='\0';

							/*Read string, then broadcast: */
							if(fread(string,string_size*sizeof(char),1,this->fid)!=1)_error_("  could not read string ");
							#ifdef _HAVE_MPI_
							MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
							#endif
						}
						else{
							string=(char*)xmalloc(sizeof(char));
							string[0]='\0';
						}
						
						/*Add string to parameters: */
						this->constants->AddObject(new StringParam(record_enum,string));

						/*Free string*/
						xfree((void**)&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_("%s%i","unknown record type:",record_code); 
						break;;
				}
			}
		}
	} //}}}
	#ifdef _HAVE_MPI_
	else{ //cpu ~0 {{{2
		for(;;){ //wait on cpu 0
			MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD);  /*get from cpu 0 what we are going to do: */
			if(record_code==0){
				break; //we are done, break from the loop
			}
			else{
				MPI_Bcast(&record_enum,1,MPI_INT,0,MPI_COMM_WORLD);   //get from cpu 0 name of the data
				MPI_Bcast(&record_length,1,MPI_INT,0,MPI_COMM_WORLD);  
				switch(record_code){
				case 1: 
					/*boolean. get it from cpu 0 */
					MPI_Bcast(&booleanint,1,MPI_INT,0,MPI_COMM_WORLD); 
						
					/*create BoolParam: */
					this->constants->AddObject(new BoolParam(record_enum,(bool)booleanint)); //cast to a boolean
					break;

				case 2:
					/*integer. get it from cpu 0 */
					MPI_Bcast(&integer,1,MPI_INT,0,MPI_COMM_WORLD); 
						
					/*create IntParam: */
					this->constants->AddObject(new IntParam(record_enum,integer));

					break;
				case 3:
					/*scalar. get it from cpu 0 */
					MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
						
					/*create DoubleParam: */
					this->constants->AddObject(new DoubleParam(record_enum,scalar));

					break;
				case 4: 
					MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
					if(string_size){
						string=(char*)xmalloc((string_size+1)*sizeof(char));
						string[string_size]='\0';

						/*Read string from cpu 0: */
						MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
					}
					else{
						string=(char*)xmalloc(sizeof(char));
						string[0]='\0';
					}
					/*Add string to parameters: */
					this->constants->AddObject(new StringParam(record_enum,string));

					/*Free string*/
					xfree((void**)&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_("%s%i","unknown record type:",record_code); 
					break;;
				}


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

	extern int my_rank;
	extern int num_procs;
	

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

	if(code!=1)_error_("expecting a boolean for enum %s",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 ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&booleanint,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

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

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

	extern int my_rank;
	extern int num_procs;
	

	/*output: */
	int   integer;
	int   code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=2)_error_("expecting an integer for enum %s",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 ");
	}

	#ifdef _HAVE_MPI_
	MPI_Bcast(&integer,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*Assign output pointers: */
	*pinteger=integer;

}
/*}}}*/
/*FUNCTION IoModel::FetchData(double*   pscalar,int data_enum){{{1*/
void  IoModel::FetchData(double* pscalar,int data_enum){


	extern int my_rank;
	extern int num_procs;
	

	/*output: */
	double   scalar;
	int      code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=3)_error_("expecting a double for enum %s",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(double),1,fid)!=1)_error_(" could not read scalar ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&scalar,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
	#endif

	/*Assign output pointers: */
	*pscalar=scalar;
		 
}
/*}}}*/
/*FUNCTION IoModel::FetchData(char**    pstring,int data_enum){{{1*/
void  IoModel::FetchData(char** pstring,int data_enum){

	extern int my_rank;
	extern int num_procs;
	

	/*output: */
	char* string=NULL;
	int   string_size;
	int code=0;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	
	if(code!=4)_error_("expecting a string for enum %s",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 ");
	}

	#ifdef _HAVE_MPI_
	MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*Now allocate string: */
	if(string_size){
		string=(char*)xmalloc((string_size+1)*sizeof(char));
		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 ");
		}
		#ifdef _HAVE_MPI_
		MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
		#endif
	}
	else{
		string=(char*)xmalloc(sizeof(char));
		string[0]='\0';
	}


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

	extern int my_rank;
	extern int num_procs;
	int i,j;

	/*output: */
	int M,N;
	double* matrix=NULL;
	int*    integer_matrix=NULL;
	int code=0;
	int vector_type=0;
	
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,&vector_type,data_enum);

	if((code!=5) && (code!=6) && (code!=7))_error_("expecting a double, integer or boolean matrix for enum %s",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 ");
	}

	#ifdef _HAVE_MPI_
	MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	if(my_rank==0){  
		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
	#endif

	/*Now allocate matrix: */
	if(M*N){
		matrix=(double*)xmalloc(M*N*sizeof(double));

		/*Read matrix on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(matrix,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix ");
		}
		
		#ifdef _HAVE_MPI_
		MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); 
		#endif
	}

	/*Now cast to integer: */
	if(M*N){
		integer_matrix=(int*)xmalloc(M*N*sizeof(int));
		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:*/
	xfree((void**)&matrix);

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

}
/*}}}*/
/*FUNCTION IoModel::FetchData(double**  pdoublematrix,int* pM,int* pN,int data_enum){{{1*/
void  IoModel::FetchData(double** pmatrix,int* pM,int* pN,int data_enum){

	extern int my_rank;
	extern int num_procs;

	/*output: */
	int M,N;
	double* matrix=NULL;
	int code=0;
	int vector_type=0;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,&vector_type,data_enum);
	if((code!=5) && (code!=6) && (code!=7))_error_("expecting a double, integer or boolean matrix for enum %s",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 ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	if(my_rank==0){  
		if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*Now allocate matrix: */
	if(M*N){
		matrix=(double*)xmalloc(M*N*sizeof(double));

		/*Read matrix on node 0, then broadcast: */
		if(my_rank==0){  
			if(fread(matrix,M*N*sizeof(double),1,fid)!=1) _error_("could not read matrix ");
		}
		#ifdef _HAVE_MPI_
		MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); 
		#endif
	}

	/*Assign output pointers: */
	*pmatrix=matrix;
	if (pM)*pM=M;
	if (pN)*pN=N;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(char***   pstrings,int* pnumstrings,int data_enum){{{1*/
void  IoModel::FetchData(char*** pstrings,int* pnumstrings,int data_enum){

	extern int my_rank;
	extern int num_procs;
	
	int i;

	/*output: */
	int   numstrings=0;
	char** strings=NULL;
	
	/*intermediary: */
	char* string=NULL;
	int   string_size;
	int   code;
	
	/*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 %s",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");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&numstrings,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	/*Now allocate string array: */
	if(numstrings){
		strings=(char**)xmalloc(numstrings*sizeof(char*));
		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 ");
			}
			#ifdef _HAVE_MPI_
			MPI_Bcast(&string_size,1,MPI_INT,0,MPI_COMM_WORLD); 
			#endif
			if(string_size){
				string=(char*)xmalloc((string_size+1)*sizeof(char));
				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 ");
				}
				#ifdef _HAVE_MPI_
				MPI_Bcast(string,string_size,MPI_CHAR,0,MPI_COMM_WORLD); 
				#endif
			}
			else{
				string=(char*)xmalloc(sizeof(char));
				string[0]='\0';
			}

			strings[i]=string;
		}
	}

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

	int i;

	extern int my_rank;
	extern int num_procs;

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

	/*intermediary: */
	int     M, N;
	double *matrix = NULL;
	int     code;
	
	/*Set file pointer to beginning of the data: */
	fid=this->SetFilePointerToData(&code,NULL,data_enum);
	if(code!=8)_error_("expecting a double mat array for enum %s",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 ");
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&numrecords,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

	if(numrecords){

		/*Allocate matrices :*/
		matrices=(double**)xmalloc(numrecords*sizeof(double*));
		mdims=(int*)xmalloc(numrecords*sizeof(int));
		ndims=(int*)xmalloc(numrecords*sizeof(int));

		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_("%s%i%s","could not read number of rows in ",i,"th matrix of matrix array");
			}
			#ifdef _HAVE_MPI_
			MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD); 
			#endif

			if(my_rank==0){  
				if(fread(&N,sizeof(int),1,fid)!=1) _error_("%s%i%s","could not read number of columns in ",i,"th matrix of matrix array");
			}
			#ifdef _HAVE_MPI_
			MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD); 
			#endif

			/*Now allocate matrix: */
			if(M*N){
				matrix=(double*)xmalloc(M*N*sizeof(double));

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

				#ifdef _HAVE_MPI_
				MPI_Bcast(matrix,M*N,MPI_DOUBLE,0,MPI_COMM_WORLD); 
				#endif
			}

			/*Assign: */
			matrices[i]=matrix;
			mdims[i]=M;
			ndims[i]=N;
		}
	}

	/*Assign output pointers: */
	*pmatrices=matrices;
	*pmdims=mdims;
	*pndims=ndims;
	*pnumrecords=numrecords;
}
/*}}}*/
/*FUNCTION IoModel::FetchData(int num,...){{{1*/
void  IoModel::FetchData(int num,...){

	va_list ap;
	int     dataenum;
	double* 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 double* : */
	
	va_start(ap,num);
	for(i=0; i<num; i++){
		
		dataenum=va_arg(ap, int);

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

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

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

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

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


	/*variables being fetched: */
	bool    boolean;
	int     integer;
	double  scalar;
	char   *string        = NULL;
	double *doublevector  = 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.  {{{2
			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=(Element*)elements->GetObjectByOffset(counter);
					element->InputCreate((double)boolean,vector_enum,code); 
					counter++;
				}
			}
			break; /*}}}*/
		case 2: //integer constant.  {{{2
			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=(Element*)elements->GetObjectByOffset(counter);
					element->InputCreate((double)integer,vector_enum,code); 
					counter++;
				}
			}
			break; /*}}}*/
		case 3: //double constant.  {{{2
			this->FetchData(&scalar,vector_enum);

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

			/*Fetch vector:*/
			this->FetchData(&doublevector,&M,&N,vector_enum); //we still have a doublevector, because it might include times in transient mode
			/*Check we got something, otherwise fetch default: */
			if(doublevector){ 
				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(&doublevector,&M,&N,default_vector_enum);
					if(doublevector){
						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=(Element*)elements->GetObjectByOffset(counter);
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(doublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}
			break; /*}}}*/
		case 6: //int vector.  {{{2

			/*Fetch vector:*/
			this->FetchData(&doublevector,&M,&N,vector_enum); //we still have a doublevector, because it might include times in transient mode
			/*Check we got something, otherwise fetch default: */
			if(doublevector){ 
				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(&doublevector,&M,&N,default_vector_enum);
					if(doublevector){
						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=(Element*)elements->GetObjectByOffset(counter);
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(doublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}
			break; /*}}}*/
		case 7: //double vector.  {{{2

			/*Fetch vector:*/
			this->FetchData(&doublevector,&M,&N,vector_enum);
			/*Check we got something, otherwise fetch default: */
			if(doublevector){ 
				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(&doublevector,&M,&N,default_vector_enum);
					if(doublevector){
						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=(Element*)elements->GetObjectByOffset(counter);
					if(defaulting) element->InputCreate(default_value,vector_enum,code); 
					else           element->InputCreate(doublevector,i,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
					counter++;
				}
			}

			break; /*}}}*/
		default: /*{{{2*/
			_error_("%s%i%s","data code ",code," not supported yet!");
			break;
			/*}}}*/
	}
	/*Free ressources:*/
	xfree((void**)&doublevector);
	xfree((void**)&string);
}
/*FUNCTION IoModel::SetFilePointerToData{{{1*/
FILE* IoModel::SetFilePointerToData(int* pcode,int* pvector_type, int data_enum){

	extern int my_rank;
	extern int num_procs;
	
	int found=0;
	int record_enum;
	int record_length;
	int record_code; //1 to 7 number
	int vector_type; //nodal or elementary
 
	/*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);
				fread(&record_code,sizeof(int),1,fid);

				/*if record_code points to a vector, get its type (nodal or elementary): */
				if(5<=record_code && record_code<=7)fread(&vector_type,sizeof(int),1,fid);
				found=1;
				break;
			}
			else{
				/*This is not the correct string, read the record length, and use it to skip this record: */
				fread(&record_length,sizeof(int),1,fid);
				/*skip: */
				fseek(fid,record_length,SEEK_CUR);
			}
		}
	}
	#ifdef _HAVE_MPI_
	MPI_Bcast(&found,1,MPI_INT,0,MPI_COMM_WORLD); 
	if(!found)_error_("%s %s ","could not find data with name",EnumToStringx(data_enum));
	#endif

	/*Broadcast code and vector type: */
	#ifdef _HAVE_MPI_
	MPI_Bcast(&record_code,1,MPI_INT,0,MPI_COMM_WORLD); 
	MPI_Bcast(&vector_type,1,MPI_INT,0,MPI_COMM_WORLD); 
	if(record_code==5) MPI_Bcast(&vector_type,1,MPI_INT,0,MPI_COMM_WORLD); 
	#endif

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

	return fid;
}
/*}}}*/
