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

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

void GroundinglineMigrationx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){

	int      i, migration_style,analysis_type;
	IssmDouble*  vertices_potentially_ungrounding = NULL;
	IssmDouble*  vertices_ungrounding             = NULL;
	IssmDouble*  old_floatingice                  = NULL;
	Vector<IssmDouble>*      vec_old_floatingice              = NULL;
	Element* element                          = NULL;

	if(VerboseModule()) _pprintLine_("   Migrating grounding line");

	/*retrieve parameters: */
	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	if(migration_style==NoneEnum) return;
	if(migration_style!=AgressiveMigrationEnum && migration_style!=SoftMigrationEnum) _error_(EnumToStringx(migration_style) << " not supported yet!");

	if(migration_style==SoftMigrationEnum){
		/*Create flag for grounded vertices above the hydrostatic equilibrium: */
		vertices_potentially_ungrounding=PotentialSheetUngrounding(elements,vertices,parameters);

		/*propagate ice shelf into connex areas of the ice sheet that potentially want to unground: */
		vertices_ungrounding=PropagateFloatingiceToGroundedNeighbors(elements,nodes,vertices,parameters,vertices_potentially_ungrounding);
	}

	/*Create vector with vertices initially floating and serialize*/
	vec_old_floatingice=CreateNodesOnFloatingIce(nodes,analysis_type);
	old_floatingice=vec_old_floatingice->ToMPISerial();

	/*Migrate grounding line : */
	for(i=0;i<elements->Size();i++){
		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
		element->MigrateGroundingLine(old_floatingice,vertices_ungrounding);
	}

	/*free ressouces: */
	xdelete(&vec_old_floatingice);
	xDelete<IssmDouble>(vertices_potentially_ungrounding);
	xDelete<IssmDouble>(vertices_ungrounding);
	xDelete<IssmDouble>(old_floatingice);
}

/*FUNCTION CreateNodesOnFloatingIce {{{*/
Vector<IssmDouble>* CreateNodesOnFloatingIce(Nodes* nodes,int configuration_type){ 

	int     i,numnods;
	Vector<IssmDouble>*   vec_nodes_on_floatingice = NULL;
	Node *node                     = NULL;

	/*First, initialize nodes_on_floatingice, which will track which nodes have changed status: */
	numnods=nodes->NumberOfNodes(configuration_type);
	vec_nodes_on_floatingice=new Vector<IssmDouble>(numnods);

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

	/*Assemble vector: */
	vec_nodes_on_floatingice->Assemble();

	return vec_nodes_on_floatingice;
}
/*%}}}*/
/*FUNCTION PotentialSheetUngrounding {{{*/
IssmDouble*    PotentialSheetUngrounding(Elements* elements,Vertices* vertices,Parameters* parameters){ 

	int      i,numberofvertices;
	IssmDouble*  vertices_potentially_ungrounding      = NULL;
	Vector<IssmDouble>*      vec_vertices_potentially_ungrounding  = NULL;
	Element* element                               = NULL;

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

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

	/*Assemble vector and serialize */
	vec_vertices_potentially_ungrounding->Assemble();
	vertices_potentially_ungrounding=vec_vertices_potentially_ungrounding->ToMPISerial();

	/*free ressouces and return: */
	xdelete(&vec_vertices_potentially_ungrounding);
	return vertices_potentially_ungrounding;
}
/*}}}*/
/*FUNCTION PropagateFloatingiceToGroundedNeighbors {{{*/
IssmDouble*    PropagateFloatingiceToGroundedNeighbors(Elements* elements,Nodes* nodes,Vertices* vertices,Parameters* parameters,IssmDouble* vertices_potentially_ungrounding){ 

	int      i,analysis_type;
	int      numberofvertices;
	int      nflipped,local_nflipped;
	IssmDouble*  nodes_on_floatingice                  = NULL;
	IssmDouble*  elements_neighboring_floatingce      = NULL;
	Vector<IssmDouble>*      vec_elements_neighboring_floatingice = NULL;
	Vector<IssmDouble>*      vec_nodes_on_floatingice              = NULL;
	Node*    node                                  = NULL;
	Element* element                               = NULL;

	/*recover parameters: */
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
	numberofvertices=vertices->NumberOfVertices();

	/*recover vec_nodes_on_floatingice*/
	vec_nodes_on_floatingice=CreateNodesOnFloatingIce(nodes,analysis_type);
	nodes_on_floatingice=vec_nodes_on_floatingice->ToMPISerial();

	nflipped=1; //bootstrap
	while(nflipped){

		/*Vector of size number of elements*/
		vec_elements_neighboring_floatingice=new Vector<IssmDouble>(elements->NumberOfElements(),true);

		/*Figure out if any of the nodes of the element will be floating -> elements neighbouting the floating ice*/
		for(i=0;i<elements->Size();i++){
			element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
			vec_elements_neighboring_floatingice->SetValue(element->Sid(),element->IsNodeOnShelfFromFlags(nodes_on_floatingice)?1.0:0.0,INS_VAL);
		}

		/*Assemble vector and serialize: */
		vec_elements_neighboring_floatingice->Assemble();
		elements_neighboring_floatingce=vec_elements_neighboring_floatingice->ToMPISerial();

		/*Go through elements_neighboring_floatingce, and update vector of the nodes that will start floating*/
		local_nflipped=0;
		for(i=0;i<elements->Size();i++){
			element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
			if(reCast<int,IssmDouble>(elements_neighboring_floatingce[element->Sid()])){
				local_nflipped+=element->UpdatePotentialSheetUngrounding(vertices_potentially_ungrounding,vec_nodes_on_floatingice,nodes_on_floatingice);
			}
		}
		vec_nodes_on_floatingice->Assemble();

		#ifdef _HAVE_MPI_
		MPI_Allreduce(&local_nflipped,&nflipped,1,MPI_INT,MPI_SUM,IssmComm::GetComm());
		if(VerboseConvergence()) _pprintLine_("   Additional number of vertices allowed to unground: " << nflipped);
		#else
		nflipped=local_nflipped;
		#endif

		/*Avoid leaks: */
		xDelete<IssmDouble>(elements_neighboring_floatingce);
		xDelete<IssmDouble>(nodes_on_floatingice);

		/*Assemble and serialize:*/
		xdelete(&vec_elements_neighboring_floatingice);
		nodes_on_floatingice=vec_nodes_on_floatingice->ToMPISerial();
	}

	/*Free ressources:*/
	xdelete(&vec_nodes_on_floatingice);
	xDelete<IssmDouble>(elements_neighboring_floatingce);

	return nodes_on_floatingice;
}
/*}}}*/
