/*!\file:  CreateEdges.cpp
 * \brief: create edges from 2d mesh
 */ 

#include "../../classes/classes.h"
#include "../../shared/shared.h"

void CreateEdges(IoModel* iomodel){

	/*If edges are already present, exit*/
	if(iomodel->edges) return;

	/*Check Iomodel properties*/
	if(iomodel->numberofvertices<3) _error_("not enough elements in mesh");
	_assert_(iomodel->elements);

	/*Intermediaries*/
	bool exist;
	int  i,j,v1,v2,v3;
	int  maxnbe,nbe,elementnbe,elementnbv;
	int *elementedges         = NULL;
	int *elementedges_markers = NULL;

	/*Mesh dependent variables*/
	if(iomodel->meshtype==Mesh2DhorizontalEnum){
		elementnbv = 3;
		elementnbe = 3;
		elementedges         = xNew<int>(elementnbe*2);
		elementedges_markers = xNew<int>(elementnbe);
		elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = 1;
		elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = 1;
		elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
	}
	else if(iomodel->meshtype==Mesh2DverticalEnum){
		elementnbv = 3;
		elementnbe = 3;
		elementedges         = xNew<int>(elementnbe*2);
		elementedges_markers = xNew<int>(elementnbe);
		elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = 1;
		elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = 1;
		elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
	}
	else if(iomodel->meshtype==Mesh3DtetrasEnum){
		elementnbv = 4;
		elementnbe = 6;
		elementedges         = xNew<int>(elementnbe*2);
		elementedges_markers = xNew<int>(elementnbe);
		elementedges[2*0+0] = 1; elementedges[2*0+1] = 2; elementedges_markers[0] = 1;
		elementedges[2*1+0] = 2; elementedges[2*1+1] = 0; elementedges_markers[1] = 1;
		elementedges[2*2+0] = 0; elementedges[2*2+1] = 1; elementedges_markers[2] = 1;
		elementedges[2*3+0] = 0; elementedges[2*3+1] = 3; elementedges_markers[3] = 1;
		elementedges[2*4+0] = 1; elementedges[2*4+1] = 3; elementedges_markers[4] = 1;
		elementedges[2*5+0] = 2; elementedges[2*5+1] = 3; elementedges_markers[5] = 1;
	}
	else if(iomodel->meshtype==Mesh3DEnum){
		elementnbv = 6;
		elementnbe = 9;
		elementedges         = xNew<int>(elementnbe*2);
		elementedges_markers = xNew<int>(elementnbe);
		elementedges[2*0+0] = 0; elementedges[2*0+1] = 3; elementedges_markers[0] = 2;
		elementedges[2*1+0] = 1; elementedges[2*1+1] = 4; elementedges_markers[1] = 2;
		elementedges[2*2+0] = 2; elementedges[2*2+1] = 5; elementedges_markers[2] = 2;
		elementedges[2*3+0] = 1; elementedges[2*3+1] = 2; elementedges_markers[3] = 1;
		elementedges[2*4+0] = 2; elementedges[2*4+1] = 0; elementedges_markers[4] = 1;
		elementedges[2*5+0] = 0; elementedges[2*5+1] = 1; elementedges_markers[5] = 1;
		elementedges[2*6+0] = 4; elementedges[2*6+1] = 5; elementedges_markers[6] = 1;
		elementedges[2*7+0] = 5; elementedges[2*7+1] = 3; elementedges_markers[7] = 1;
		elementedges[2*8+0] = 3; elementedges[2*8+1] = 4; elementedges_markers[8] = 1;
	}
	else{
		_error_("mesh dimension not supported yet");
	}

	/*Maximum number of edges*/
	maxnbe = elementnbe*iomodel->numberofelements;

	/*Initialize intermediaries*/
	int *edgestemp                 = xNew<int>(maxnbe*3);                             /*format: [vertex1 vertex2 marker]       */
	int *element_edge_connectivity = xNew<int>(iomodel->numberofelements*elementnbe); /*format: [edge1 edge2 ... edgen] */

	/*Initialize chain*/
	int* head_minv = xNew<int>(iomodel->numberofvertices);
	int* next_edge = xNew<int>(maxnbe);
	for(i=0;i<iomodel->numberofvertices;i++) head_minv[i]=-1;

	/*Initialize number of edges*/
	nbe = 0;

	for(i=0;i<iomodel->numberofelements;i++){
		for(j=0;j<elementnbe;j++){

			/*Get the two indices of the edge number j of the ith element*/
			v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
			v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);

			/*v1 and v2 must be sorted*/
			if(v2<v1){
				v3=v2; v2=v1; v1=v3;
			}

			/*This edge a priori has not been processed yet*/
			exist = false;

			/*Go through all processed edges connected to v1 and check whether we have seen this edge yet*/
			for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
				if(edgestemp[e*3+1]==v2+1){
					exist = true;
					element_edge_connectivity[i*elementnbe+j]=e;
					break;
				}
			}

			/*If this edge is new, add it to the lists*/
			if(!exist){
				_assert_(nbe<maxnbe);

				/*Update edges*/
				edgestemp[nbe*3+0] = v1+1;
				edgestemp[nbe*3+1] = v2+1;
				edgestemp[nbe*3+2] = elementedges_markers[j];

				/*Update Connectivity*/
				element_edge_connectivity[i*elementnbe+j]=nbe;

				/*Update chain*/
				next_edge[nbe] = head_minv[v1];
				head_minv[v1]  = nbe;

				/*Increase number of edges*/
				nbe++;
			}
		}
	}

	/*Clean up*/
	xDelete<int>(head_minv);
	xDelete<int>(next_edge);
	xDelete<int>(elementedges_markers);

	/*Create final edges*/
	int* edges = xNew<int>(nbe*3); /*format: [vertex1 vertex2 marker]*/
	for(int i=0;i<3*nbe;i++) edges[i] = edgestemp[i];

	/*Clean up*/
	xDelete<int>(edgestemp);
	xDelete<int>(elementedges);

	/*Assign output pointers*/
	iomodel->edges                     = edges;
	iomodel->elementtoedgeconnectivity = element_edge_connectivity;
	iomodel->numberofedges             = nbe;
}
