/*
 * CreateConstraintsEnthalpy.c:
 */

#include "../../../Container/Container.h"
#include "../../../io/io.h"
#include "../../../toolkits/toolkits.h"
#include "../../../EnumDefinitions/EnumDefinitions.h"
#include "../../../objects/objects.h"
#include "../../../shared/shared.h"
#include "../../../include/include.h"
#include "../ModelProcessorx.h"

void	CreateConstraintsEnthalpy(Constraints** pconstraints, IoModel* iomodel){

	/*Intermediary*/
	int    i,j;
	int    count;
	int    dim;
	int    M,N;
	int    numberofvertices;
	bool   spcpresent=false;
	double heatcapacity;
	double referencetemperature;
	
	/*Output*/
	IssmDouble *spcvector  = NULL;
	IssmDouble* times=NULL;
	IssmDouble* values=NULL;
	Constraints* constraints = NULL;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
	iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);

	/*Recover pointer: */
	constraints=*pconstraints;

	/*Create constraints if they do not exist yet*/
	if(!constraints) constraints = new Constraints();

	/*return if 2d mesh*/
	if (dim==2){
		*pconstraints=constraints;
		return;
	}

	/*Fetch data: */
	iomodel->FetchData(&spcvector,&M,&N,ThermalSpctemperatureEnum);

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

		for (i=0;i<numberofvertices;i++){
			/*keep only this partition's nodes:*/
			if((iomodel->my_vertices[i])){

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

					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,heatcapacity*(spcvector[i]-referencetemperature),EnthalpyAnalysisEnum));
					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]=spcvector[(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]=heatcapacity*(spcvector[i*N+j]-referencetemperature);
					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,EnthalpyAnalysisEnum));
					count++;
				}
				xDelete<IssmDouble>(values);
			}
		}
	}
	else{
		_error2_("Size of field " << EnumToStringx(ThermalSpctemperatureEnum) << " not supported");
	}

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

	/*Assign output pointer: */
	*pconstraints=constraints;
}
