/*!\file GroundinglineMigrationx
 * \brief: migration grounding line position.
 */

#include "./GroundinglineMigrationx.h"
#include "./GroundinglineMigrationxLocal.h"
#include "../../shared/shared.h"
#include "../../io/io.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"

/*FUNCTION PotentialSheetUngrounding {{{1*/
double*    PotentialSheetUngrounding(Elements* elements,Vertices* vertices,Parameters* parameters){ 

	int      i,numberofvertices;
	double*  vertices_potentially_ungrounding      = NULL;
	Vec      vec_vertices_potentially_ungrounding  = NULL;
	Element* element                               = NULL;

	/*Initialize vector with number of vertices*/
	numberofvertices=vertices->NumberOfVertices();
	vec_vertices_potentially_ungrounding=NewVec(numberofvertices); //grounded vertex that could start floating

	/*Fill vector vertices_potentially_floating: */
	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		element->PotentialSheetUngrounding(vec_vertices_potentially_ungrounding);
	}

	/*Assemble vector: */
	VecAssemblyBegin(vec_vertices_potentially_ungrounding);
	VecAssemblyEnd(vec_vertices_potentially_ungrounding);

	/*Serialize vector: */
	VecToMPISerial(&vertices_potentially_ungrounding,vec_vertices_potentially_ungrounding);

	/*free ressouces and return: */
	VecFree(&vec_vertices_potentially_ungrounding);
	return vertices_potentially_ungrounding;
}
/*}}}*/
/*FUNCTION PropagateShelfIntoConnexIceSheet {{{1*/
double*    PropagateShelfIntoConnexIceSheet(Elements* elements,Nodes* nodes,Parameters* parameters,double* potential_sheet_ungrounding){ 

	int      i;
	Element* element=NULL;
	int      numnods;
	int      analysis_type;
	Vec      vec_nodes_on_iceshelf=NULL;
	double*  nodes_on_iceshelf=NULL;
	int      nflipped,local_nflipped;
	double*  elements_touching_iceshelf=NULL;

	/*recover parameters: */
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*recover vec_nodes_on_iceshelf: */
	vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,analysis_type);

	nflipped=1; //bootstrap
	while(nflipped){
		
		/*get a list of potential elements that have nodes on ice shelf: */
		elements_touching_iceshelf=CreateElementTouchingIceShelf(elements,vec_nodes_on_iceshelf);

		/*now,  go through elements_touching_iceshelf, and if they have nodes inside potential_sheet_ungrounding, 
		 * flag it: */
		local_nflipped=0;

		/*serialize vec_nodes_on_iceshelf, needed by element's UpdatePotentialSheetUngrounding routine, to figure out if 
		 * nodes have flipped from grounded to ungrounded: */
		VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);

		for(i=0;i<elements->Size();i++){
			element=(Element*)elements->GetObjectByOffset(i);
			if(elements_touching_iceshelf[element->Sid()]){
				local_nflipped+=element->UpdatePotentialSheetUngrounding(potential_sheet_ungrounding,vec_nodes_on_iceshelf,nodes_on_iceshelf);
			}
		}
		
		MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
		_printf_(VerboseConvergence(),"   number of grounded vertices  connected to grounding line: %i\n",nflipped);

		/*Avoid leaks: */
		xfree((void**)&nodes_on_iceshelf);
		xfree((void**)&elements_touching_iceshelf);

		/*Assemble:*/
		VecAssemblyBegin(vec_nodes_on_iceshelf);
		VecAssemblyEnd(vec_nodes_on_iceshelf);
	}

	/*Serialize: */
	VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);

	/*Free ressources:*/
	VecFree(&vec_nodes_on_iceshelf);
	xfree((void**)&elements_touching_iceshelf);

	return nodes_on_iceshelf;
}
/*}}}*/
/*FUNCTION CreateElementOnGroundingline {{{1*/
bool*      CreateElementOnGroundingline(Elements* elements,double* element_on_iceshelf){ 

	int      i,j;
	Element *element       = NULL;
	bool    *element_on_gl = NULL;
	bool     ongl=false;
	int     *neighboorsids = NULL;
	int      sid;
	int      shelf;

	/*Go through elements, and look for elements that can possibly have a grounding line migration. These 
	 * are  elements that touch the grounding line: */
	element_on_gl=(bool*)xmalloc(elements->Size()*sizeof(bool));

	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		
		sid=element->Sid();
		neighboorsids=element->GetHorizontalNeighboorSids();

		shelf=(int)element_on_iceshelf[sid];
		ongl=false;
		for(j=0;j<3;j++){
			if (neighboorsids[j]<0)continue; //this neighboor does not exist
			if ((shelf==1) & (element_on_iceshelf[neighboorsids[j]]==1))continue;  //both this element and this neighboor are on th ice shelf
			if ((shelf==0) & (element_on_iceshelf[neighboorsids[j]]==0))continue;  //both this element and this neighboor are on the ice sheet
			ongl=true; //neighboor j is on a different location than us, ie we are touching the grounding line.
		}
		element_on_gl[i]=ongl;
	}

	return element_on_gl;
}
/*}}}*/
/*FUNCTION CreateElementOnIceShelf {{{1*/
double*    CreateElementOnIceShelf(Elements* elements){ 

	int i;
	Element* element=NULL;
	Vec vec_element_on_iceshelf=NULL;
	double* element_on_iceshelf=NULL;

	/*Create  vector holding  all the elements IsFloating flags: */
	vec_element_on_iceshelf=NewVec(elements->NumberOfElements(),true);

	/*Loop through elements, and fill vec_element_on_iceshelf: */
	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		VecSetValue(vec_element_on_iceshelf,element->Sid(),(int)element->IsFloating(),INSERT_VALUES);
	}

	/*Assemble vector: */
	VecAssemblyBegin(vec_element_on_iceshelf);
	VecAssemblyEnd(vec_element_on_iceshelf);

	/*Serialize vector: */
	VecToMPISerial(&element_on_iceshelf,vec_element_on_iceshelf);

	/*free ressouces: */
	VecFree(&vec_element_on_iceshelf);

	return element_on_iceshelf;
}
/*}}}*/
/*FUNCTION CreateElementTouchingIceShelf {{{1*/
double*    CreateElementTouchingIceShelf(Elements* elements,Vec vec_nodes_on_iceshelf){ 

	int i;
	Element* element=NULL;
	Vec vec_element_touching_iceshelf=NULL;
	double* nodes_on_iceshelf=NULL;
	
	/*output: */
	double* element_touching_iceshelf=NULL;

	/*serialize: */
	VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);

	/*Create  vector holding  all the elements IsFloating flags: */
	vec_element_touching_iceshelf=NewVec(elements->NumberOfElements(),true);

	/*Loop through elements, and fill vec_element_touching_iceshelf: */
	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		VecSetValue(vec_element_touching_iceshelf,element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_iceshelf)?1.0:0.0,INSERT_VALUES);
	}

	/*Assemble vector: */
	VecAssemblyBegin(vec_element_touching_iceshelf);
	VecAssemblyEnd(vec_element_touching_iceshelf);

	/*Serialize vector: */
	VecToMPISerial(&element_touching_iceshelf,vec_element_touching_iceshelf);

	/*free ressouces: */
	VecFree(&vec_element_touching_iceshelf);
	xfree((void**)&nodes_on_iceshelf);

	return element_touching_iceshelf;
}
/*}}}*/
/*FUNCTION UpdateShelfStatus {{{1*/
int        UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){

	int i;
	Element* element=NULL;
	int configuration_type;
	int numnods;
	Vec vec_new_shelf_nodes=NULL;
	double* new_shelf_nodes=NULL;

	/*output: */
	int local_nflipped=0;
	int nflipped=0;

	/*recover parameters: */
	parameters->FindParam(&configuration_type,ConfigurationTypeEnum);

	/*First, initialize vec_new_shelf_nodes, which will track which nodes have changed status: */
	numnods=nodes->NumberOfNodes(configuration_type);
	vec_new_shelf_nodes=NewVec(numnods);

	/*Ok, now go through  the elements that have gone through grounding line migration, 
	 * and update their flags: */
	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		if(element_touching_iceshelf[element->Sid()]){
			local_nflipped+=element->UpdateShelfStatus(vec_new_shelf_nodes);
		}
	}
	MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);

	/*Serialize vec_new_shelf_nodes: */
	VecToMPISerial(&new_shelf_nodes,vec_new_shelf_nodes);

	/*Now, go through ALL elements, and update the status of the nodes, to propagate what happened at the grounding line: */
	/*Carry out grounding line migration for those elements: */
	for(i=0;i<elements->Size();i++){
		element=(Element*)elements->GetObjectByOffset(i);
		element->UpdateShelfFlags(new_shelf_nodes);
	}

	/*Free ressources: */
	VecFree(&vec_new_shelf_nodes);
	xfree((void**)&new_shelf_nodes);

	return nflipped;
}
/*}}}*/
/*FUNCTION CreateNodesOnIceShelf {{{1*/
Vec CreateNodesOnIceShelf(Nodes* nodes,int configuration_type){ 

	/*output: */
	Vec     nodes_on_iceshelf=NULL;

	/*intermediary: */
	int     numnods;
	int     i;
	Node*   node=NULL;

	/*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
	numnods=nodes->NumberOfNodes(configuration_type);
	nodes_on_iceshelf=NewVec(numnods);

	/*Loop through nodes, and fill nodes_on_iceshelf: */
	for(i=0;i<nodes->Size();i++){
		node=(Node*)nodes->GetObjectByOffset(i);
		if(node->InAnalysis(configuration_type)){
			if(node->IsFloating()){
				VecSetValue(nodes_on_iceshelf,node->Sid(),1.0,INSERT_VALUES);
			}
		}
	}

	/*Assemble vector: */
	VecAssemblyBegin(nodes_on_iceshelf);
	VecAssemblyEnd(nodes_on_iceshelf);

	return nodes_on_iceshelf;
}
/*}}}*/