/*! \file CreateLoadsStressbalance.c:
 */
#include <string.h>
#include "../../../toolkits/toolkits.h"
#include "../../../classes/classes.h"
#include "../../../shared/shared.h"
#include "../ModelProcessorx.h"

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

	/*DataSets*/
	Pengrid   *pengrid   = NULL;

	/*Intermediary*/
	int         segment_width;
	int         element;
	int         i;
	int         count;
	int         penpair_ids[2];
	bool        isSSA,isL1L2,isHO,isFS;
	int         numpenalties,numberofpressureloads,numrifts,numriftsegments;
	IssmDouble *pressureload   = NULL;
	IssmDouble *elements_type  = NULL;
	IssmDouble *nodeoniceshelf = NULL;
	IssmDouble *riftinfo       = NULL;
	IssmDouble *nodeonbed      = NULL;
	IssmDouble *nodeonFS   = NULL;
	IssmDouble *nodeonicesheet = NULL;
	IssmDouble *vertices_type  = NULL;
	IssmDouble *penalties      = NULL;
	int         assert_int;

	/*Fetch parameters: */
	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
	iomodel->Constant(&isFS,FlowequationIsFSEnum);
	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
	iomodel->Constant(&isHO,FlowequationIsHOEnum);
	iomodel->Constant(&numrifts,RiftsNumriftsEnum);

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

	/*Now, is the flag macayaealHO on? otherwise, do nothing: */
	if(!isSSA & !isHO & !isFS & !isL1L2) return;

	/*Initialize counter: */
	count=0;

	/*Create Penpair for penalties: */
	iomodel->FetchData(&penalties,&numpenalties,NULL,StressbalanceVertexPairingEnum);

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

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

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

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

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

	/*free ressources: */
	iomodel->DeleteData(penalties,StressbalanceVertexPairingEnum);

	/*Create Riffront loads for rifts: */
	#ifdef _HAVE_RIFTS_
	if(numrifts){
		iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum);
		iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
		for(i=0;i<numriftsegments;i++){
			if(iomodel->my_elements[reCast<int,IssmDouble>(*(riftinfo+RIFTINFOSIZE*i+2))-1]){
				loads->AddObject(new Riftfront(iomodel->loadcounter+count+1,i,iomodel,StressbalanceAnalysisEnum));
				count++;
			}
		}
		iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskGroundediceLevelsetEnum);
		xDelete<IssmDouble>(riftinfo);
	}
	#endif

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