/*!\file IoModelToConstraintsx
 * \brief: used in ModelProcessor, to retrieve a nodal vector  and to create corresponding constraints using Constraints objects
 */

#include "./IoModelToConstraintsx.h"
#include "../../shared/shared.h"
#include "../../io/io.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type){

	/*intermediary: */
	int     i,j;
	bool        transient        = false;
	FILE       *fid              = NULL;
	int         code             = 0;
	int         vector_layout    = 0;
	IssmDouble *times            = NULL;
	IssmDouble *values           = NULL;
	bool        spcpresent       = false;
	int         count            = 0;
	int         numberofvertices;

	/*variables being fetched: */
	IssmDouble *IssmDoublevector  = NULL;
	int     M,N;

	/*Fetch parameters: */
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);

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

	if(code!=7)_error_("expecting a IssmDouble vector for constraints with enum " << EnumToStringx(vector_enum));
	if(vector_layout!=1)_error_("expecting a nodal vector for constraints with enum " << EnumToStringx(vector_enum));

	/*Fetch vector:*/
	iomodel->FetchData(&IssmDoublevector,&M,&N,vector_enum);

	/*Transient or static?:*/
	if(M==numberofvertices){
		/*static: just create Constraints objects*/
		count=0;

		/*Create Constraints from x,y,z: */
		for (i=0;i<numberofvertices;i++){

			/*keep only this partition's nodes:*/
			if((iomodel->my_vertices[i])){

				if (!xIsNan<IssmDouble>(IssmDoublevector[i])){

					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,IssmDoublevector[i],analysis_type));
					count++;
				}
			}
		}
	}
	else if (M==(numberofvertices+1)){
		/*transient: create transient SpcTransient objects. Same logic, except we need to retrieve 
		 * various times and values to initialize an SpcTransient object: */
		count=0;

		/*figure out times: */
		times=xNew<IssmDouble>(N);
		for(j=0;j<N;j++){
			times[j]=IssmDoublevector[(M-1)*N+j];
		}
		/*unit conversion: */
		UnitConversion(times,N,ExtToIuEnum,TimeEnum);

		/*Create constraints from x,y,z: */
		for (i=0;i<numberofvertices;i++){

			/*keep only this partition's nodes:*/
			if((iomodel->my_vertices[i])){

				/*figure out times and values: */
				values=xNew<IssmDouble>(N);
				spcpresent=false;
				for(j=0;j<N;j++){
					values[j]=IssmDoublevector[i*N+j];
					if(!xIsNan<IssmDouble>(values[j]))spcpresent=true; //NaN means no spc by default
				}

				if(spcpresent){
					constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,N,times,values,analysis_type));
					count++;
				}
				xDelete<IssmDouble>(values);
			}
		}
	}
	else{
		_error_("Size of field " << EnumToStringx(vector_enum) << " not supported");
	}

	/*Free ressources:*/
	xDelete<IssmDouble>(IssmDoublevector);
	xDelete<IssmDouble>(times);
	xDelete<IssmDouble>(values);
}
