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

#undef __FUNCT__ 
#define __FUNCT__ "SplitMeshForRifts"

#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 grid;
	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: */
	RiftSegmentsFromSegments(&nriftsegs,&riftsegments,nel,index,nsegs,segments); /*riftsegments of size nriftsegsx4 (4 for first element on segment,second element, 
																				   first grid and second sgrid)*/

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

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

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

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

	cleanup_and_return:

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