/*!\file:  CreateSingleNodeToElementConnectivity.cpp
 * \brief: create connectivity table
 */ 

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

#include "../../shared/shared.h"
#include "../../io/io.h"
#include "../../include/include.h"
#include "./ModelProcessorx.h"

void CreateSingleNodeToElementConnectivity(IoModel* iomodel){

	/*Intermediary*/
	int i,j;
	int vertexid;
	int elementswidth;
	int    dim;
	int    numberofelements;
	int    numberofvertices;
	IssmDouble* elements=NULL;

	/*output*/
	int* connectivity=NULL;

	/*Return if connectivity already present*/
	if(iomodel->singlenodetoelementconnectivity) return;

	/*Fetch parameters: */
	iomodel->Constant(&dim,MeshDimensionEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	elements=iomodel->Data(MeshElementsEnum);

	/*Some checks if debugging*/
	_assert_(numberofvertices);
	_assert_(numberofelements);
	_assert_(iomodel->my_elements);
	_assert_(elements);

	/*Allocate ouput*/
	connectivity=xNewZeroInit<int>(numberofvertices);

	/*Get element width (3 or 6)*/
	if (dim==2){
		elementswidth=3;
	}
	else{
		elementswidth=6;
	}

	/*Create connectivity table*/
	for (i=0;i<numberofelements;i++){
		/*!! in parallel we do not want the vertex to be connected to an element that is not in its partition!!*/
		if(iomodel->my_elements[i]){
			for (j=0;j<elementswidth;j++){
				vertexid=reCast<int>(elements[elementswidth*i+j]);
				_assert_(vertexid>0 && vertexid-1<numberofvertices);
				connectivity[vertexid-1]=i+1;
			}
		}
	}

	/*Assign to iomodel*/
	iomodel->singlenodetoelementconnectivity=connectivity;

	/*Do not free connectivity!!!*/
}
