/*! \file CreateLoadsDiagnosticHoriz.c:
 */

#undef __FUNCT__ 
#define __FUNCT__ "CreateLoadsDiagnosticHoriz"

#include "../DataSet/DataSet.h"
#include "../toolkits/toolkits.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../objects/objects.h"
#include "../shared/shared.h"
#include "../include/macros.h"
#include "./ModelProcessorx.h"


int	CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){

	int noerr=1;
	int i,j,counter;
	int element;

	extern int my_rank;
	extern int num_procs;
	
	DataSet*    loads    = NULL;
	Icefront*   icefront = NULL;
	Penpair*    penpair  = NULL;

	int segment_width;
	int element_type;
	int i1,i2,i3,i4;
	int i0,range;
	int grid1,grid2;

	/*rift penpair: */
	double  normal[2];
	double  length;
	double  friction;
	int     fill;
		
	/*icefront intermediary data: */
	char icefront_type[ICEFRONTSTRING];
	int icefront_element_type;
	int	icefront_sid;
	int icefront_eid;
	int icefront_mparid;
	int	icefront_node_ids[MAX_ICEFRONT_GRIDS];
	double icefront_h[MAX_ICEFRONT_GRIDS];
	double	icefront_b[MAX_ICEFRONT_GRIDS];

	/*penpair intermediary data: */
	int penpair_id;
	int penpair_node_ids[2];
	double penpair_penalty_offset;
	int penpair_numdofs;
	int penpair_dof;
	int penpair_penalty_lock;
	int penpair_element_ids[2];
	double penpair_friction;
	int penpair_fill;
	double penpair_normal[2];
	double penpair_length;

	/*Rifts:*/
	int* riftsnumpenaltypairs=NULL;
	double** riftspenaltypairs=NULL;
	int* riftsfill=NULL;
	int* riftsfriction=NULL;
	double* riftpenaltypairs=NULL;
	int el1,el2;
	
	/*Create loads: */
	loads   = new DataSet(LoadsEnum());
	
	/*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
	 * referenced by a certain load must belong to the cluster node): */
	ModelFetchData((void**)&model->segmentonneumann_diag,&model->numberofsegs_diag,NULL,model_handle,"segmentonneumann_diag","Matrix","Mat");
	ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");

	/*First load data:*/
	for (i=0;i<model->numberofsegs_diag;i++){
		
		if (strcmp(model->meshtype,"2d")==0){
			segment_width=3;
			element_type=TriaEnum();
		}
		else{
			segment_width=5;
			element_type=PentaEnum();
		}


		element=(int)(*(model->segmentonneumann_diag+segment_width*i+segment_width-1)-1); //element is in the last column

		#ifdef _PARALLEL_
		if (model->epart[element]!=my_rank){
			/*This load does not belong to this cluster node, as it references an element 
			 *that does not belong to this node's partition. Drop this 'i':*/
			continue;
		}
		#endif
	
		icefront_mparid=model->numberofelements+1; //matlab indexing
		icefront_sid=i+1; //matlab indexing
		icefront_eid=(int)*(model->segmentonneumann_diag+segment_width*i+segment_width-1); //matlab indexing
		icefront_element_type=element_type;

		i1=(int)*(model->segmentonneumann_diag+segment_width*i+0);
		i2=(int)*(model->segmentonneumann_diag+segment_width*i+1);
			
		icefront_node_ids[0]=i1;
		icefront_node_ids[1]=i2;
		
		icefront_h[0]=model->thickness[i1-1];
		icefront_h[1]=model->thickness[i2-1];

		icefront_b[0]=model->bed[i1-1];
		icefront_b[1]=model->bed[i2-1];
		
		if ((int)*(model->elements_type+2*element+0)==MacAyealEnum()){ //this is a collapsed 3d element, icefront will be 2d
			strcpy(icefront_type,"segment");
		}
		else if ((int)*(model->elements_type+2*element+0)==PattynEnum()){ //this is a real 3d element, icefront will be 3d.
			strcpy(icefront_type,"quad");
			i3=(int)*(model->segmentonneumann_diag+segment_width*i+2);
			i4=(int)*(model->segmentonneumann_diag+segment_width*i+3);
			icefront_node_ids[2]=i3;
			icefront_node_ids[3]=i4;
			
			icefront_h[2]=model->thickness[i3-1];
			icefront_h[3]=model->thickness[i4-1];

			icefront_b[2]=model->bed[i3-1];
			icefront_b[3]=model->bed[i4-1];
		}
		else{
			_printf_("%s%s\n",__FUNCT__," error message: element type not supported yet");
			return 0;
		}

		icefront = new Icefront(icefront_type,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b);
		
		loads->AddObject(icefront);

	}
	/*Free data: */
	xfree((void**)&model->segmentonneumann_diag);
	xfree((void**)&model->elements_type);
	xfree((void**)&model->thickness);
	xfree((void**)&model->bed);

		
	/*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */

	/*First fetch penalties: */
	if (strcmp(model->meshtype,"3d")==0){
		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
	}

	counter=0;
	if(model->numpenalties){

		/*First deal with internal grids: */
		if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU

			for(j=0;j<model->numlayers-1;j++){

				/*We are pairing grids along a vertical profile.*/
				grid1=(int)*(model->penalties+model->numlayers*i+j);
				grid2=(int)*(model->penalties+model->numlayers*i+j+1);
				
				penpair_id=counter+1; //matlab indexing
				penpair_dof=1;
				penpair_node_ids[0]=grid1;
				penpair_node_ids[1]=grid2;
				penpair_penalty_offset=model->penalty_offset;
				penpair_numdofs=1;
		
				penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
						penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
				loads->AddObject(penpair);
		

				counter++;
			
				penpair_id=counter+1; //matlab indexing
				penpair_dof=2;
				penpair_node_ids[0]=grid1;
				penpair_node_ids[1]=grid2;
				penpair_penalty_offset=model->penalty_offset;
				penpair_numdofs=1;
			
				penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
						penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
				loads->AddObject(penpair);
	

				counter++;

			}
		} //for ( i=0; i<numpenalties; i++) 
	}

	/*Free data: */
	xfree((void**)&model->penalties);
	
	/*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
	/*First fetch rifts: */
	if(model->numrifts){

		for(i=0;i<model->numrifts;i++){
			riftpenaltypairs=model->riftspenaltypairs[i];
			for(j=0;j<model->riftsnumpenaltypairs[i];j++){

				el1=(int)*(riftpenaltypairs+7*j+2);
				#ifdef _PARALLEL_
				if (model->epart[el1-1]!=my_rank){
					/*This load does not belong to this cluster node, as it references an element 
					 *that does not belong to this node's partition. Drop this 'j':*/
					continue;
				}
				#endif

				/*Ok, retrieve all the data needed to add a penalty between the two grids: */
				el2=(int)*(riftpenaltypairs+7*j+3); 
				grid1=(int)*(riftpenaltypairs+7*j+0); 
				grid2=(int)*(riftpenaltypairs+7*j+1);
				normal[0]=*(riftpenaltypairs+7*j+4);
				normal[1]=*(riftpenaltypairs+7*j+5);
				length=*(riftpenaltypairs+7*j+6);
				friction=model->riftsfriction[i];
				fill=model->riftsfill[i];

				penpair_id=counter+1; //matlab indexing
				penpair_node_ids[0]=grid1;
				penpair_node_ids[1]=grid2;
				penpair_element_ids[0]=el1;
				penpair_element_ids[1]=el2;
				penpair_penalty_offset=model->penalty_offset;
				penpair_penalty_lock=model->penalty_lock;
				penpair_numdofs=2;
				penpair_normal[0]=normal[0];
				penpair_normal[1]=normal[1];
				penpair_length=length;
				penpair_friction=friction;
				penpair_fill=fill;

				penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
						penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
				
				loads->AddObject(penpair);

				counter++;
			}
		}
	}

	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
	 * datasets, it will not be redone: */
	loads->Presort();

	/*Free ressources:*/
	xfree((void**)&riftsnumpenaltypairs); 
	for(i=0;i<model->numrifts;i++){
		double* temp=riftspenaltypairs[i];
		xfree((void**)&temp);
	}
	xfree((void**)&riftspenaltypairs);
	xfree((void**)&riftsfill);
	xfree((void**)&riftsfriction);


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

	return noerr;
}


