/*! \file CreateLoadsPrognostic.c:
 */

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

void	CreateLoadsPrognostic(Loads** ploads, IoModel* iomodel){

	/*Intermediaries*/
	int i;
	int element;
	int penpair_ids[2];
	int count=0;
	int stabilization;
	int numberofedges;
	int numvertex_pairing;
	int numberofelements;

	/*DataSet*/
	Loads*    loads    = NULL;
	
	/*Fetch parameters: */
	iomodel->Constant(&stabilization,PrognosticStabilizationEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&numberofedges,MeshNumberofedgesEnum);

	/*Recover pointer: */
	loads=*ploads;

	/*Create loads if they do not exist yet*/
	if(!loads) loads = new Loads();

	/*Loads only in DG*/
	if (stabilization==3){

		/*Get edges and elements*/
		iomodel->FetchData(3,MeshEdgesEnum,MeshElementsEnum,ThicknessEnum);

		/*First load data:*/
		for (i=0;i<numberofedges;i++){

			/*Get left and right elements*/
			element=reCast<int>(iomodel->Data(MeshEdgesEnum)[4*i+2])-1; //edges are [node1 node2 elem1 elem2]

			/*Now, if this element is not in the partition, pass: */
			if(!iomodel->my_elements[element]) continue;

			/* Add load */
			loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,iomodel,PrognosticAnalysisEnum));
		}

		/*Free data: */
		iomodel->DeleteData(3,MeshEdgesEnum,MeshElementsEnum,ThicknessEnum);
	}

	/*Create Penpair for vertex_pairing: */
	IssmDouble *vertex_pairing=NULL;
	IssmDouble *nodeonbed=NULL;
	iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,PrognosticVertexPairingEnum);
	iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);

	for(i=0;i<numvertex_pairing;i++){

		if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){

			/*In debugging mode, check that the second node is in the same cpu*/
			_assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);

			/*Skip if one of the two is not on the bed*/
			if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;

			/*Get node ids*/
			penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
			penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);

			/*Create Load*/
			loads->AddObject(new Penpair(
							iomodel->loadcounter+count+1,
							&penpair_ids[0],
							PrognosticAnalysisEnum));
			count++;
		}
	}

	/*free ressources: */
	xDelete<IssmDouble>(vertex_pairing);
	xDelete<IssmDouble>(nodeonbed);

	/*Assign output pointer: */
	*ploads=loads;
}
