/*
 * SplitMeshForRifts.c:
 */
#include "./trimesh.h"

#include "../Alloc/alloc.h"

int SplitMeshForRifts(int* pnel,double** pindex,int* pnods,double** px,double** py,int* pnsegs,double** psegments,double** psegmentmarkerlist){

	/*Some notes on dimensions: 
	index  of size nelx3
	x and y of size nodsx1
	segments of size nsegsx3*/
	
	/*Error management: */
	int noerr=1;
	
	int i,j,k,l;
	int node;
	int el;

	int  nriftsegs;
	int* riftsegments=NULL; 
	int* flags=NULL;
	
	int  NumGridElementListOnOneSideOfRift;
	int* GridElementListOnOneSideOfRift=NULL;

	/*Input: */
	int     nel;
	double* index=NULL;
	int     nods;
	double* x=NULL;
	double* y=NULL;
	double* segments=NULL;
	double*    segmentmarkerlist=NULL;
	int     nsegs;

	/*Recover input: */
	nel=*pnel;
	index=*pindex;
	nods=*pnods;
	x=*px;
	y=*py;
	nsegs=*pnsegs;
	segments=*psegments;
	segmentmarkerlist=*psegmentmarkerlist;

	/*Establish list of segments that belong to a rift: */
	/*riftsegments of size nriftsegsx4 (4 for first element on segment,second element,first node and second snode)*/
	RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments);

	/*Go through all nodes of the rift segments, and start splitting the mesh: */
	flags=(int*)xcalloc(nods,sizeof(int)); //to make sure we don't split the same nodes twice!
	for (i=0;i<nriftsegs;i++){
		for (j=0;j<2;j++){
	
			node=riftsegments[4*i+j+2];
			if(flags[node-1]){
				/*This node was already split, skip:*/
				continue;
			}
			else{
				flags[node-1]=1;
			}

			if(IsGridOnRift(riftsegments,nriftsegs,node)){
			
				DetermineGridElementListOnOneSideOfRift(&NumGridElementListOnOneSideOfRift,&GridElementListOnOneSideOfRift,i,nriftsegs,riftsegments,node,index,nel);
			
				/*Summary: we have for node, a list of elements
				 * (GridElementListOnOneSideOfRift, of size
				 * NumGridElementListOnOneSideOfRift) that all contain node 
				 *and that are on the same side of the rift. For all these
				 elements, we clone node into another node, and we swap all
				 instances of node in the triangulation *for those elements, to the
				 new node.*/
				
				//augment number of nodes 
				nods=nods+1;
				//create new node
				x=(double*)xrealloc(x,nods*sizeof(double));
				y=(double*)xrealloc(y,nods*sizeof(double));
				x[nods-1]=x[node-1]; //matlab indexing
				y[nods-1]=y[node-1]; //matlab indexing

				//change elements owning this node
				for (k=0;k<NumGridElementListOnOneSideOfRift;k++){
					el=GridElementListOnOneSideOfRift[k];
					for (l=0;l<3;l++){
						if (*(index+3*el+l)==node)*(index+3*el+l)=nods; //again, matlab indexing.
					}
				}
			}// if(IsGridOnRift(riftsegments,nriftsegs,node))
		} //for(j=0;j<2;j++)
	} //for (i=0;i<nriftsegs;i++)

	/*update segments: they got modified completely by adding new nodes.*/
	UpdateSegments(&segments,&segmentmarkerlist, &nsegs,index,x,y,riftsegments,nriftsegs,nods,nel);

	/*Assign output pointers: */
	*pnel=nel;
	*pindex=index;
	*pnods=nods;
	*px=x;
	*py=y;
	*pnsegs=nsegs;
	*psegments=segments;
	*psegmentmarkerlist=segmentmarkerlist;
	return noerr;
}
