/*!\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 CreateElementOnGroundingline {{{1*/
bool*      CreateElementOnGroundingline(Elements* elements,double* element_on_iceshelf){ 

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

	/*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;
	double  *element_on_iceshelf     = NULL;
	Element *element                 = NULL;
	Vec      vec_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,double* nodes_on_iceshelf){ 

	int      i;
	double  *element_touching_iceshelf     = NULL;
	Element *element                       = NULL;
	Vec      vec_element_touching_iceshelf = NULL;

	/*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);

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

	int     i,numberofvertices;
	Vec     nodes_on_iceshelf = NULL;
	Node*   node              = NULL;

	/*First, initialize nodes_on_iceshelf, which will track which nodes have changed status: */
	numberofvertices=vertices->NumberOfVertices();
	nodes_on_iceshelf=NewVec(numberofvertices);

	/*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;
}
/*}}}*/
/*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,Vertices* vertices,Parameters* parameters,double* vertices_potentially_ungrounding){ 

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

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

	/*recover vec_nodes_on_iceshelf and serialize: */
	vec_nodes_on_iceshelf=CreateNodesOnIceShelf(nodes,vertices,analysis_type);
	VecToMPISerial(&nodes_on_iceshelf,vec_nodes_on_iceshelf);

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

		/*Go through elements_touching_iceshelf, and flag their nodes that can unground inside potential_sheet_ungrounding*/
		local_nflipped=0;
		for(i=0;i<elements->Size();i++){
			element=(Element*)elements->GetObjectByOffset(i);
			if(elements_touching_iceshelf[element->Sid()]){
				local_nflipped+=element->UpdatePotentialSheetUngrounding(vertices_potentially_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**)&elements_touching_iceshelf);
		xfree((void**)&nodes_on_iceshelf);

		/*Assemble and serialize:*/
		VecAssemblyBegin(vec_nodes_on_iceshelf);
		VecAssemblyEnd(vec_nodes_on_iceshelf);
		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 UpdateShelfStatus {{{1*/
int UpdateShelfStatus(Elements* elements,Nodes* nodes,Parameters* parameters,double* element_touching_iceshelf){

	int      i,numnods,configuration_type;
	int      nflipped=0;
	int      local_nflipped=0;
	Vec      vec_new_shelf_nodes = NULL;
	double*  new_shelf_nodes     = NULL;
	Element* element             = NULL;

	/*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;
}
/*}}}*/