/*!\file:  ElementsAndVerticesPartitioning.cpp
 * \brief: partition elements and nodes and vertices
 */ 

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <string.h>
#include "../../objects/objects.h"
#include "../../shared/shared.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../io/io.h"
#include "../../include/include.h"
#include "../MeshPartitionx/MeshPartitionx.h"
#include "../ModelProcessorx/ModelProcessorx.h"

void  ElementsAndVerticesPartitioning(bool** pmy_elements, int** pmy_vertices, IoModel* iomodel){

	int i;

	extern int my_rank;
	extern int num_procs;
	int    numberofelements;
	int    numberofvertices;
	int    numberofelements2d;
	int    numberofvertices2d;
	int    numlayers;
	int    numrifts;
	int    numvertex_pairing;

	/*output: */
	bool* my_elements=NULL;
	int* my_vertices=NULL;

	/*intermediary: */
	int* epart=NULL; //element partitioning.
	int* npart=NULL; //node partitioning.
	int  elements_width; //number of columns in elements (2d->3, 3d->6)
	int  el1,el2;
	int    dim;
	IssmDouble* elements=NULL;
	IssmDouble* elements2d=NULL;
	IssmDouble* riftinfo=NULL;
	IssmDouble* vertex_pairing=NULL;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	iomodel->Constant(&numberofelements2d,MeshNumberofelements2dEnum);
	iomodel->Constant(&numberofvertices2d,MeshNumberofvertices2dEnum);
	iomodel->Constant(&numlayers,MeshNumberoflayersEnum);
	iomodel->Constant(&numrifts,RiftsNumriftsEnum);

	/*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */
	if (*pmy_elements)return;

	/*Number of vertices per elements, needed to correctly retrieve data: */
	if(dim==2) elements_width=3; //tria elements
	else elements_width=6; //penta elements

	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
	if(dim==2){
		/*load elements: */
		iomodel->FetchData(&elements,NULL,NULL,MeshElementsEnum);
	}
	else{
		/*load elements2d: */
		iomodel->FetchData(&elements2d,NULL,NULL,MeshElements2dEnum);
	}

	MeshPartitionx(&epart, &npart,numberofelements,numberofvertices,elements, numberofelements2d,numberofvertices2d,elements2d,numlayers,elements_width, dim,num_procs);

	/*Free elements and elements2d: */
	xDelete<IssmDouble>(elements);
	xDelete<IssmDouble>(elements2d);

	/*Deal with rifts, they have to be included into one partition only, not several: */
	if(numrifts){
		iomodel->FetchData(&riftinfo,&numrifts,NULL,RiftsRiftstructEnum);
		for(i=0;i<numrifts;i++){
			el1=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+2))-1; //matlab indexing to c indexing
			el2=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+3))-1; //matlab indexing to c indexing
			epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding vertices;
		}
		xDelete<IssmDouble>(riftinfo); 
	}

	/*Used later on: */
	my_vertices=xNewZeroInit<int>(numberofvertices);
	my_elements=xNewZeroInit<bool>(numberofelements);

	/*Start figuring out, out of the partition, which elements belong to this cpu: */
	iomodel->FetchData(&elements,NULL,NULL,MeshElementsEnum);
	for (i=0;i<numberofelements;i++){

		/*!All elements have been partitioned above, only deal with elements for this cpu: */
		if(my_rank==epart[i]){ 

			my_elements[i]=true;
			
			/*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 
			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
			 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
			 will hold which vertices belong to this partition*/
			my_vertices[reCast<int>(*(elements+elements_width*i+0))-1]=1;
			my_vertices[reCast<int>(*(elements+elements_width*i+1))-1]=1;
			my_vertices[reCast<int>(*(elements+elements_width*i+2))-1]=1;
			
			if(elements_width==6){
				my_vertices[reCast<int>(*(elements+elements_width*i+3))-1]=1;
				my_vertices[reCast<int>(*(elements+elements_width*i+4))-1]=1;
				my_vertices[reCast<int>(*(elements+elements_width*i+5))-1]=1;
			}
		}
	}//for (i=0;i<numberofelements;i++)
	/*Free data : */
	xDelete<IssmDouble>(elements);

	/*We might have vertex_pairing in which case, some vertices have to be cloned:
	 * penpair has 2 nodes that are poointing toward 2 vertices.
	 * The 2 vertices must be in the same cpu as the penpair*/
	iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,DiagnosticVertexPairingEnum);
	for(i=0;i<numvertex_pairing;i++){
		if(my_vertices[reCast<int>(vertex_pairing[2*i+0])-1] && !my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]){
			my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]=2; //to know that these elements are not on the partition
		}
	}
	xDelete<IssmDouble>(vertex_pairing);
	iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,PrognosticVertexPairingEnum);
	for(i=0;i<numvertex_pairing;i++){
		if(my_vertices[reCast<int>(vertex_pairing[2*i+0])-1] && !my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]){
			my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]=2; //to know that these elements are not on the partition
		}
	}
	xDelete<IssmDouble>(vertex_pairing);

	/*Free ressources:*/
	xDelete<int>(npart);
	xDelete<int>(epart);

	/*Assign output pointers:*/
	*pmy_elements=my_elements;
	*pmy_vertices=my_vertices;
}
