/*! \file CreateLoadsDiagnosticHoriz.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	CreateLoadsDiagnosticHoriz(Loads** ploads, IoModel* iomodel){

	/*DataSets*/
	Icefront  *icefront  = NULL;
	Pengrid   *pengrid   = NULL;

	/*Intermediary*/
	int segment_width;
	int element;
	int i;
	int count=0;
	int penpair_ids[2];
	int dim;
	int numberofvertices;
	bool ismacayealpattyn,isstokes,isl1l2;
	int  numpenalties,numberofpressureloads,numrifts,numriftsegments;
	IssmDouble *pressureload   = NULL;
	IssmDouble *elements_type  = NULL;
	IssmDouble *nodeoniceshelf = NULL;
	IssmDouble *riftinfo       = NULL;
	IssmDouble *nodeonbed      = NULL;
	IssmDouble *nodeonstokes   = NULL;
	IssmDouble *nodeonicesheet = NULL;
	IssmDouble *vertices_type  = NULL;
	IssmDouble *penalties      = NULL;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum);
	iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
	iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
	iomodel->Constant(&numrifts,RiftsNumriftsEnum);

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

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

	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
	if(!ismacayealpattyn & !isstokes & !isl1l2){
		*ploads=loads;
		return;
	}

	/*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the nodes
	 * referenced by a certain load must belong to the cluster node): */
	iomodel->FetchData(3,DiagnosticIcefrontEnum,ThicknessEnum,BedEnum);
	iomodel->FetchData(&pressureload,&numberofpressureloads,NULL,DiagnosticIcefrontEnum);
	iomodel->FetchData(&elements_type,NULL,NULL,FlowequationElementEquationEnum);

	/*Initialize counter: */
	count=0;

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

		/*Retrieve element to which this icefront belongs: */
		if (dim==2) segment_width=4; 
		else segment_width=6;
		element=reCast<int,IssmDouble>(*(pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)

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

		/*Do not create ice front if Hutter or Stokes elements*/
		if (reCast<int,IssmDouble>(*(elements_type+element))==HutterApproximationEnum) continue;

		/*Create and  add load: */
		if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealApproximationEnum) && dim==2){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealApproximationEnum) && dim==3){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(L1L2ApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(PattynApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(L1L2ApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,L1L2IceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(StokesApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealPattynApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(PattynStokesApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if (reCast<int,IssmDouble>(*(elements_type+element))==(MacAyealStokesApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
	}

	/*Free data: */
	iomodel->DeleteData(3,DiagnosticIcefrontEnum,ThicknessEnum,BedEnum);
	iomodel->DeleteData(elements_type,FlowequationElementEquationEnum);
	iomodel->DeleteData(pressureload,DiagnosticIcefrontEnum);

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

	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_(iomodel->my_vertices[(int)penalties[2*i+1]-1]);

			/*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],DiagnosticHorizAnalysisEnum));
			count++;
		}
	}

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

	/*Create Riffront loads for rifts: */
	if(numrifts){
		iomodel->FetchData(&riftinfo,&numriftsegments,NULL,RiftsRiftstructEnum);
		iomodel->FetchData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskVertexonfloatingiceEnum);
		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,DiagnosticHorizAnalysisEnum));
				count++;
			}
		}
		iomodel->DeleteData(5,RiftsRiftstructEnum,ThicknessEnum,BedEnum,SurfaceEnum,MaskVertexonfloatingiceEnum);
		xDelete<IssmDouble>(riftinfo);
	}

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