/*! \file CreateLoadsDiagnosticHoriz.c:
 */

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

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

	/*DataSets*/
	Loads     *loads     = NULL;
	Icefront  *icefront  = NULL;
	Riftfront *riftfront = 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;
	int  numpenalties,numberofpressureloads,numrifts,numriftsegments;
	double *pressureload   = NULL;
	double *elements_type  = NULL;
	double *nodeoniceshelf = NULL;
	double *riftinfo       = NULL;
	double *nodeonbed      = NULL;
	double *nodeonstokes   = NULL;
	double *nodeonicesheet = NULL;
	double *vertices_type  = NULL;
	double *penalties      = NULL;

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

	/*Recover pointer: */
	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){
		*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=(int)(*(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 ((int)*(elements_type+element)==HutterApproximationEnum) continue;

		/*Create and  add load: */
		if ((int)*(elements_type+element)==(MacAyealApproximationEnum) && dim==2){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if ((int)*(elements_type+element)==(MacAyealApproximationEnum) && dim==3){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if ((int)*(elements_type+element)==(PattynApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if ((int)*(elements_type+element)==(StokesApproximationEnum)){
			loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum));
			count++;
		}
		else if ((int)*(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 ((int)*(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 ((int)*(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);
	xDelete<double>(elements_type);
	xDelete<double>(pressureload);

	/*Create Penpair for penalties: */
	iomodel->FetchData(&penalties,&numpenalties,NULL,DiagnosticVertexPairingEnum);
	
	for(i=0;i<numpenalties;i++){

		if(iomodel->my_vertices[(int)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+(int)penalties[2*i+0];
			penpair_ids[1]=iomodel->nodecounter+(int)penalties[2*i+1];

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

	/*free ressources: */
	xDelete<double>(penalties);

	/*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[(int)*(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<Riftfront>(riftfront);
	}

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