/*
 * TriMeshUtils: mesh manipulation routines: 
 */

#include "./trimesh.h"

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


int IsGridOnRift(int* riftsegments, int nriftsegs, int grid){

	/*Does this grid belong to 4 elements, or just 2? If it belongs to 4 elements, it is inside a rift, 
	 *if it belongs to 2 elements, it is on the tip of a rift, or it has already been split across the rift (see below).*/
	
	int i;
	int j;
	int count;

	count=0;
	for (i=0;i<nriftsegs;i++){
		for (j=0;j<2;j++){
			if ((*(riftsegments+4*i+2+j))==grid) count++;
		}
	}
	if (count==2){
		return 1;
	}
	else{
		return 0;
	}
}
				

int GridElementsList(int** pGridElements, int* pNumGridElements,int grid,double * index,int nel){

	/*From a grid, recover all the elements that are connected to it: */
	int i,j;
	int noerr=1;

	int max_number_elements=12;
	int current_size;
	int NumGridElements;
	int* GridElements=NULL;
	int* GridElementsRealloc=NULL;

	/*From a mesh with 30 degrees minimum angle, we get 12 possible elements that own 
	 * the grid. We start by allocating GridElements with that size, and realloc 
	 * more if needed.*/

	current_size=max_number_elements;
	NumGridElements=0;
	GridElements=(int*)xmalloc(max_number_elements*sizeof(int));

	for (i=0;i<nel;i++){
		for (j=0;j<3;j++){
			if ((int)*(index+3*i+j)==grid){
				if (NumGridElements<=(current_size-1)){
					GridElements[NumGridElements]=i;
					NumGridElements++;
					break;
				}
				else{
					/*Reallocate another max_number_elements slots in the GridElements: */
					GridElementsRealloc=(int*)xrealloc(GridElements,(current_size+max_number_elements)*sizeof(int));
					if (!GridElementsRealloc){
						noerr=0;
						goto cleanup_and_return;
					}
					current_size+=max_number_elements;
					GridElements=GridElementsRealloc;
					GridElements[NumGridElements]=i;
					NumGridElements++;
					break;
				}
			}
		}
	}
	cleanup_and_return:
	if(!noerr){
		xfree((void**)&GridElements);
	}
	/*Allocate return pointers: */
	*pGridElements=GridElements;
	*pNumGridElements=NumGridElements;
	return noerr;
}


int IsNeighbor(int el1,int el2,double* index){
	/*From a triangulation held in index, figure out if elements 1 and 2 have two grids in common: */
	int i,j;
	int count=0;
	for (i=0;i<3;i++){
		for (j=0;j<3;j++){
			if (*(index+3*el1+i)==*(index+3*el2+j))count++;
		}
	}
	if (count==2){
		return 1;
	}
	else{
		return 0;
	}
}
							

int IsOnRift(int el,int nriftsegs,int* riftsegments){
	/*From a list of elements segments, figure out if el belongs to it: */
	int i;
	for (i=0;i<nriftsegs;i++){
		if ((*(riftsegments+4*i+0)==el) | (*(riftsegments+4*i+1)==el)){
			return 1;
		}
	}
	return 0;
}


/******************************************************************************************************************************
                                   RiftSegmentsFromSegments
******************************************************************************************************************************/

int RiftSegmentsFromSegments(int* pnriftsegs, int** priftsegments, int nel, double* index, int nsegs,double* segments){
	
	int i,counter;
	int el,el2;
	
	int  nriftsegs;
	int* riftsegments=NULL;
	int* riftsegments_uncompressed=NULL; 
	double element_grids[3];

	/*Allocate segmentflags: */
	riftsegments_uncompressed=(int*)xcalloc(nsegs*5,sizeof(int));

	/*Find the segments that belong to a rift: they are the ones that see two elements. The other ones belong to a boundary 
	 *or a hole: */
	nriftsegs=0;
	for (i=0;i<nsegs;i++){
		el=(int)*(segments+3*i+2)-1; //element found in AssociateSegmentToElements
		/*Temporarily set grids belonging to the segments to -1 in the triangulation index, and 
		 *then  proceed to find another element that owns the segment. If we don't find it, we know 
		 *we are dealing with a boundary or hole, otherwise, we are dealing with a rift: */
		element_grids[0]=*(index+3*el+0);
		element_grids[1]=*(index+3*el+1);
		element_grids[2]=*(index+3*el+2);

		*(index+3*el+0)=-1;
		*(index+3*el+1)=-1;
		*(index+3*el+2)=-1;

		el2=FindElement(*(segments+3*i+0),*(segments+3*i+1),index,nel); 

		/*Restore index: */
		*(index+3*el+0)=element_grids[0];
		*(index+3*el+1)=element_grids[1];
		*(index+3*el+2)=element_grids[2];

		if (el2!=-1){
			/*el and el2 are on a segment rift, facing one another, plug them into riftsegments_uncompressed: */
			#ifdef _DEBUG_
				printf("Elements %i and %i are on a rift\n",el+1,el2+1);
			#endif
		    *(riftsegments_uncompressed+5*i+0)=1;
		    *(riftsegments_uncompressed+5*i+1)=el;
		    *(riftsegments_uncompressed+5*i+2)=el2;
		    *(riftsegments_uncompressed+5*i+3)=(int)*(segments+3*i+0);
		    *(riftsegments_uncompressed+5*i+4)=(int)*(segments+3*i+1);
			nriftsegs++;
		}
	}

	/*Compress riftsegments_uncompressed:*/
	riftsegments=(int*)xmalloc(nriftsegs*4*sizeof(int));
	counter=0;
	for (i=0;i<nsegs;i++){
		if (*(riftsegments_uncompressed+5*i+0)){
			*(riftsegments+counter*4+0)=*(riftsegments_uncompressed+5*i+1);
			*(riftsegments+counter*4+1)=*(riftsegments_uncompressed+5*i+2);
			*(riftsegments+counter*4+2)=*(riftsegments_uncompressed+5*i+3);
			*(riftsegments+counter*4+3)=*(riftsegments_uncompressed+5*i+4);
			counter++;
		}
	}

 	#ifdef _DEBUG_
	for (i=0;i<nriftsegs;i++){
		printf("Elements %i and %i are on a rift, joined by grids %i and %i\n",*(riftsegments+4*i+0)+1,*(riftsegments+4*i+1)+1,*(riftsegments+4*i+2),*(riftsegments+4*i+3));
	}
	#endif
	xfree((void**)&riftsegments_uncompressed);
	
	/*Assign output pointers: */
	*priftsegments=riftsegments;
	*pnriftsegs=nriftsegs;
}

/******************************************************************************************************************************
                                   DetermineGridElementListOnOneSideOfRift
******************************************************************************************************************************/

int DetermineGridElementListOnOneSideOfRift(int* pNumGridElementListOnOneSideOfRift, int** pGridElementListOnOneSideOfRift, int segmentnumber, int nriftsegs, int* riftsegments, int grid,double* index,int nel){

	int noerr=1;
	int k,l,counter;
	int newel;

	int* GridElements=NULL;
	int  NumGridElements;

	/*Output: */
	int NumGridElementListOnOneSideOfRift;
	int* GridElementListOnOneSideOfRift=NULL;

	/*Build a list of all the elements connected to this grid: */
	GridElementsList(&GridElements,&NumGridElements,grid,index,nel);
	#ifdef _DEBUG_
		printf("Connected elements for grid %i\n",grid);
		for (k=0;k<NumGridElements;k++){
			printf("El %i\n",GridElements[k]+1);
		}
	#endif
			
	/*Figure out the list of elements  that are on the same side of the rift. To do so, we start from one 
	 * side of the rift and keep rotating in the same direction:*/
	GridElementListOnOneSideOfRift=(int*)xmalloc(NumGridElements*sizeof(int));
	//bootstrap the GridElementListOnOneSideOfRift by filling elements from riftsegments: */
	GridElementListOnOneSideOfRift[0]=*(riftsegments+4*segmentnumber+0); /*this one does not belong to the same side, but is just there 
															   for a rotation direction, we 'll take it out when we are 
															   done rotating*/
	GridElementListOnOneSideOfRift[1]=*(riftsegments+4*segmentnumber+1);
	#ifdef _DEBUG_
	printf("Starting with elements %i and %i for grid %i\n",GridElementListOnOneSideOfRift[0]+1,GridElementListOnOneSideOfRift[1]+1,grid);
	#endif
	counter=1;
	for (;;){
		/*Find neighbour of element GridElementListOnOneSideOfRift[counter], not 
		 * equal to GridElementListOnOneSideOfRift[counter-1]*/
		for (k=0;k<NumGridElements;k++){
			#ifdef _DEBUG_
			printf("k: %i GridElements[k]: %i GridElementListOnOneSideOfRift[counter-1]: %i GridElementListOnOneSideOfRift[counter] %i Neighboor %i\n",k,GridElements[k]+1,
				GridElementListOnOneSideOfRift[counter-1]+1,GridElementListOnOneSideOfRift[counter]+1,IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index));
			#endif
			if(IsNeighbor(GridElements[k],GridElementListOnOneSideOfRift[counter],index)){
				/*Verify this element is not already in our list of element on the same side of the rift: */
				newel=1;
				for (l=0;l<=counter;l++){
					if (GridElements[k]==GridElementListOnOneSideOfRift[l]){
						newel=0;
						break;
					}
				}
				if (newel){
					counter++;
					GridElementListOnOneSideOfRift[counter]=GridElements[k];
					if (IsOnRift(GridElements[k],nriftsegs,riftsegments)){
						break;
					}
					k=-1;
				}
			}
		}
		/*Reduce counter by 1 and get rift of first element in GridElementListOnOneSideOfRift:*/
		NumGridElementListOnOneSideOfRift=counter;
		for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
			GridElementListOnOneSideOfRift[l]=GridElementListOnOneSideOfRift[l+1];
		}
		#ifdef _DEBUG_
		printf("Grid %i is owned by  the following elements  on the same side of the rift: \n",grid);
		for (l=0;l<NumGridElementListOnOneSideOfRift;l++){
			printf("%i ",GridElementListOnOneSideOfRift[l]+1);
		}
		printf("\n");
		#endif
		break;
	}// for (;;)

	/*Free ressources: */
	xfree((void**)&GridElements);
	/*Assign output pointers: */
	*pNumGridElementListOnOneSideOfRift=NumGridElementListOnOneSideOfRift;
	*pGridElementListOnOneSideOfRift=GridElementListOnOneSideOfRift;
	return noerr;
}

/******************************************************************************************************************************
                                   UpdateSegments
******************************************************************************************************************************/

int UpdateSegments(double** psegments,double** psegmentmarkerlist, int* pnsegs, double* index, double* x,double* y,int* riftsegments,int nriftsegs){

	int noerr=1;
	int i,j,k;
	int el1,el2;

	double* segments=NULL;
	double*    segmentmarkerlist=NULL;
	int     nsegs;

	/*Recover input: */
	segments=*psegments;
	segmentmarkerlist=*psegmentmarkerlist;
	nsegs=*pnsegs;

	/*Reallocate segments: */
	segments=(double*)xrealloc(segments,(nsegs+nriftsegs)*3*sizeof(double));
	segmentmarkerlist=(double*)xrealloc(segmentmarkerlist,(nsegs+nriftsegs)*sizeof(double));

	/*First, update the existing segments to the new grids :*/
	for (i=0;i<nriftsegs;i++){
		el1=*(riftsegments+4*i+0);
		el2=*(riftsegments+4*i+1);
		for (j=0;j<nsegs;j++){
			if (*(segments+3*j+2)==(el1+1)){
				#ifdef _DEBUG_
					printf("Segment %i is the same as rift segment %i\n",j,i);
				#endif
				/*segment j is the same as rift segment i.Let's update segments[j][:] using  element el1 and the corresponding rift segment.
				 *Because riftsegments does not represent a list of rift segments anymore (it got heavily modified in SplitElementsForRifts, 
				 *we can only rely on the position (x,y) of the rift grids to create a segment:*/
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+0)-1]) & (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
						*(segments+3*j+0)=*(index+el1*3+k);
						break;
					}
				}
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+1)-1]) & (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
						*(segments+3*j+1)=*(index+el1*3+k);
						break;
					}
				}
				/*Deal with el2: */
				*(segments+3*(nsegs+i)+2)=el2+1;
				*(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+0)-1]) & (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
						*(segments+3*(nsegs+i)+0)=*(index+el2*3+k);
						break;
					}
				}
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+1)-1]) & (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
						*(segments+3*(nsegs+i)+1)=*(index+el2*3+k);
						break;
					}
				}
			}
			if (*(segments+3*j+2)==(el2+1)){
				#ifdef _DEBUG_
					printf("Segment %i is the same as rift segment %i\n",j,i);
				#endif
				/*segment j is the same as rift segment i.*/
				/*Let's update segments[j][:] using  element el2 and the corresponding rift segment: */
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+0)-1]) & (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
						*(segments+3*j+0)=*(index+el2*3+k);
						break;
					}
				}
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el2*3+k)-1]==x[(int)*(segments+3*j+1)-1]) & (y[(int)*(index+el2*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
						*(segments+3*j+1)=*(index+el2*3+k);
						break;
					}
				}
				/*Deal with el1: */
				*(segments+3*(nsegs+i)+2)=el1+1;
				*(segmentmarkerlist+(nsegs+i))=*(segmentmarkerlist+j);
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+0)-1]) & (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+0)-1])){
						*(segments+3*(nsegs+i)+0)=*(index+el1*3+k);
						break;
					}
				}
				for (k=0;k<3;k++){
					if ((x[(int)*(index+el1*3+k)-1]==x[(int)*(segments+3*j+1)-1]) & (y[(int)*(index+el1*3+k)-1]==y[(int)*(segments+3*j+1)-1])){
						*(segments+3*(nsegs+i)+1)=*(index+el1*3+k);
						break;
					}
				}
			}
		}
	}
	nsegs+=nriftsegs;

	/*Assign output pointers: */
	*psegments=segments;
	*psegmentmarkerlist=segmentmarkerlist;
	*pnsegs=nsegs;
	
	return noerr;
}

/******************************************************************************************************************************
                                   pnpoly
******************************************************************************************************************************/
int pnpoly(int npol, double *xp, double *yp, double x, double y) {
	int i, j, c = 0;
	for (i = 0, j = npol-1; i < npol; j = i++) {
		if ((((yp[i]<=y) && (y<yp[j])) ||
					((yp[j]<=y) && (y<yp[i]))) &&
				(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
			c = !c;
	}
	return c;
}

/******************************************************************************************************************************
                                   IsInPoly
******************************************************************************************************************************/
int IsInPoly(double* in,double* xc,double* yc,int numgrids,double* x,double* y,int nods){

	int i,j;
	double x0,y0;

	/*Go through all grids of the mesh:*/
	for (i=0;i<nods;i++){
		if (in[i]){
			/*this grid already is inside one of the contours, continue*/
			continue;
		}
		/*pick up grid: */
		x0=x[i];
		y0=y[i];
		if (pnpoly(numgrids,xc,yc,x0,y0)){
			in[i]=1;
		}
	}
}

/******************************************************************************************************************************
                                   FindElement
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "FindElement"

int FindElement(double A,double B,double* index,int nel){

	int n;
	int el=-1;
	#ifdef _DEBUG_
		printf("Looking for %lf %lf\n",A,B);
	#endif
	for (n=0;n<nel;n++){
		if (((*(index+3*n+0)==A)  || (*(index+3*n+1)==A) || (*(index+3*n+2)==A) ) && ((*(index+3*n+0)==B)  || (*(index+3*n+1)==B) || (*(index+3*n+2)==B))){
			el=n;
			#ifdef _DEBUG_
				printf("Found them: %lf %lf %lf\n",*(index+3*n+0),*(index+3*n+1),*(index+3*n+2));
			#endif
			break;
		}
	}
	return el;
}
/******************************************************************************************************************************
                                   SplitRiftSegments
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "SplitRiftSegments"

int SplitRiftSegments(double** psegments,double** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,double*** priftssegments,int numrifts){

	/*Using segment markers, wring out the rift segments from the segments. Rift markers are 
	 *of the form 2+i where i=0 to number of rifts */

	int noerr=1;
	int i,j,counter;

	/*input: */
	double* segments=NULL;
	double* segmentmarkerlist=NULL;
	int numsegs;
	
	/*output: */
	int* riftsnumsegs=NULL;
	double** riftssegments=NULL;
	int new_numsegs;
	double* new_segments=NULL;
	double* new_segmentmarkers=NULL;

	/*intermediary: */
	double* riftsegment=NULL;

	/*Recover input arguments: */
	segments=*psegments;
	numsegs=*pnumsegs;
	segmentmarkerlist=*psegmentmarkerlist;

	/*First, figure out  how many segments will be left in 'segments': */
	counter=0;
	for (i=0;i<numsegs;i++){
		if (segmentmarkerlist[i]==1)counter++; //1 is default marker for non-rifts;
	}
	/*Allocate new segments: */
	new_numsegs=counter;
	new_segments=(double*)xmalloc(new_numsegs*3*sizeof(double));
	new_segmentmarkers=(double*)xmalloc(new_numsegs*sizeof(double));

	/*Copy new segments info : */
	counter=0;
	for (i=0;i<numsegs;i++){
		if (segmentmarkerlist[i]==1){
			*(new_segments+3*counter+0)=*(segments+3*i+0);
			*(new_segments+3*counter+1)=*(segments+3*i+1);
			*(new_segments+3*counter+2)=*(segments+3*i+2);
			new_segmentmarkers[counter]=segmentmarkerlist[i];
			counter++;
		}
	}

	/*Now deal with rift segments: */
	riftsnumsegs=(int*)xmalloc(numrifts*sizeof(int));
	riftssegments=(double**)xmalloc(numrifts*sizeof(double*));
	for (i=0;i<numrifts;i++){
		/*Figure out how many segments for rift i: */
		counter=0;
		for (j=0;j<numsegs;j++){
			if (segmentmarkerlist[j]==2+i)counter++;
		}
		riftsnumsegs[i]=counter;
		riftsegment=(double*)xmalloc(counter*3*sizeof(double));
		/*Copy new segments info :*/
		counter=0;
		for (j=0;j<numsegs;j++){
			if (segmentmarkerlist[j]==(2+i)){
				*(riftsegment+3*counter+0)=*(segments+3*j+0);
				*(riftsegment+3*counter+1)=*(segments+3*j+1);
				*(riftsegment+3*counter+2)=*(segments+3*j+2);
				counter++;
			}
		}
		*(riftssegments+i)=riftsegment;
	}

	/*Free ressources: */
	xfree((void**)&segments);

	/*Assign output pointers: */
	*psegments=new_segments;
	*psegmentmarkerlist=new_segmentmarkers;
	*pnumsegs=new_numsegs;
	*pnumrifts=numrifts;
	*priftssegments=riftssegments;
	*priftsnumsegs=riftsnumsegs;
	return noerr;
}

/******************************************************************************************************************************
                                   PairRiftElements
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "PairRiftElements"
		
		
int PairRiftElements(int** priftsnumpairs, double*** priftspairs,int numrifts,int* riftsnumsegments, double** riftssegments,double* x,double* y){


	int noerr=1;
	int i,j,k;

	/*output: */
	int* riftsnumpairs=NULL;
	double** riftspairs=NULL;

	/*intermediary :*/
	int     numsegs;
	double* segments=NULL;
	double* pairs=NULL;
	int     grid1,grid2,grid3,grid4;

	riftsnumpairs=(int*)xmalloc(numrifts*sizeof(int));
	riftspairs=(double**)xmalloc(numrifts*sizeof(double*));
	for (i=0;i<numrifts;i++){
		segments=riftssegments[i];
		numsegs=riftsnumsegments[i];
		riftsnumpairs[i]=numsegs;
		pairs=(double*)xmalloc(2*numsegs*sizeof(double));
		for (j=0;j<numsegs;j++){
			*(pairs+2*j+0)=*(segments+3*j+2); //retrieve element to which this segment belongs.
			grid1=(int)*(segments+3*j+0)-1; grid2=(int)*(segments+3*j+1)-1;
			/*Find element facing on other side of rift: */
			for (k=0;k<numsegs;k++){
				if (k==j)continue;
				grid3=(int)*(segments+3*k+0)-1; grid4=(int)*(segments+3*k+1)-1;
				/*We are trying to find 2 elements, where position of grid3 == position of grid1, and position of grid4 == position of grid2*/
				if (   (x[grid3]==x[grid1]) && (y[grid3]==y[grid1]) && (x[grid4]==x[grid2]) && (y[grid4]==y[grid2])
				    || (x[grid3]==x[grid2]) && (y[grid3]==y[grid2]) && (x[grid4]==x[grid1]) && (y[grid4]==y[grid1])  ){
					/*We found the corresponding element: */
					*(pairs+2*j+1)=*(segments+3*k+2);
					break;
				}
			}
		}
		riftspairs[i]=pairs;
	}


	/*Assign output pointers: */
	*priftsnumpairs=riftsnumpairs;
	*priftspairs=riftspairs;

	return noerr;
}


/******************************************************************************************************************************
                                   RemoveRifts
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "RemoveRifts"
		
double dabs(double x){
	if (x<0)x=-x;
	return x;
}
int RemoveRifts(double** pindex,double** px,double** py,int* pnods,double** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,double** rifts1segments,double** rifts1pairs,int nel){

	int noerr=1;
	int i,j,k,counter,counter1,counter2;

	/*input: */
	double* index=NULL;
	double* x=NULL;
	double* y=NULL;
	int     nods;
	double* segments=NULL;
	int     numsegs;

	/*intermediary: */
	double* riftsegments=NULL;
	double* riftpairs=NULL;
	int     grid1,grid2,grid3,grid4,temp_grid;
	double  el1,el2;
	int     newnods; //temporary # grid counter.
	double  xmin,ymin;
	double* xreal=NULL;
	double* yreal=NULL;
	int* grids=NULL;
	int* merginggrids=NULL;
	int     max_size;
	int     redundant;


	/*Recover input: */
	index=*pindex;
	x=*px;
	y=*py;
	nods=*pnods;;
	segments=*psegments;
	numsegs=*pnumsegs;

	/*initialize newnods : */
	newnods=nods;

	/*Figure out a unique value to flag x and y for grid removal: */
	xmin=x[0];
	ymin=y[0];
	for (i=0;i<nods;i++){
		if (x[i]<xmin)xmin=x[i];
		if (y[i]<ymin)ymin=y[i];
	}
	xmin=xmin-dabs(xmin); 
	ymin=ymin-dabs(ymin);

	/*Initialize two arrays, one for grids that are going to be merged, the other with corresponding grids being merge into: */
	max_size=0;
	for (i=0;i<numrifts1;i++){
		max_size+=rifts1numsegs[i];
	}
	grids=(int*)xmalloc(max_size*sizeof(int));
	merginggrids=(int*)xmalloc(max_size*sizeof(int));

	/*Go through the rifts segments, and identify which grid we are going to merge with its counterpart on the other side 
	 *of the rift. The way we identify this grid is by looking at the element pairs, and the corresponding grids: */
	counter=0;
	for (i=0;i<numrifts1;i++){
		riftsegments=rifts1segments[i];
		riftpairs=rifts1pairs[i];
		for (j=0;j<rifts1numsegs[i];j++){
			el1=*(riftpairs+2*j+0);
			el2=*(riftpairs+2*j+1);
			grid1=(int)*(riftsegments+3*j+0);
			grid2=(int)*(riftsegments+3*j+1);
			/*Summary, el1 and el2 are facing one another across the rift. grid1 and grid2 belong to el1 and 
			 *are located on the rift. Find grid3 and grid4, grids belonging to el2 and located on the rift: */
			for (k=0;k<rifts1numsegs[i];k++){
				if (*(riftsegments+3*k+2)==el2){
					grid3=(int)*(riftsegments+3*k+0);
					grid4=(int)*(riftsegments+3*k+1);
					break;
				}
			}
			/* Make sure grid3 faces grid1 and grid4 faces grid2: */
			if ((x[grid1]==x[grid4]) && (y[grid1]==y[grid4])){
				/*Swap grid3 and grid4:*/
				temp_grid=grid3;
				grid3=grid4;
				grid4=temp_grid;
			}
			/* Is any of these two grid pairs on the tip of a rift, in which case, we don't include it in grids or merginggrids: */
			if ((grid1==grid3) || (grid2==grid4)){
				if(grid1!=grid3){
					/*Add grid1 and grid3 to grids and merginggrids if they have not already been added: */
					redundant=0;
					for (k=0;k<counter;k++){
						if ((merginggrids[k]==grid1) || (grids[k]==grid1))redundant=1;
					}
					if(!redundant){
						/*Ok, add grid1 to grids, and grid3 to merginggrids: */
						grids[counter]=grid1;
						merginggrids[counter]=grid3;
						counter++;
					}
				}
				if(grid2!=grid4){
					/*Add grid2 and grid4 to grids and merginggrids if they have not already been added: */
					redundant=0;
					for (k=0;k<counter;k++){
						if ((merginggrids[k]==grid2) || (grids[k]==grid2))redundant=1;
					}
					if(!redundant){
						/*Ok, add grid2 to grids, and grid4 to merginggrids: */
						grids[counter]=grid2;
						merginggrids[counter]=grid4;
						counter++;
					}
				}
			}
			else{
				/*Check that grid1 is not already present in the merginggrids: */
				redundant=0;
				for (k=0;k<counter;k++){
					if ((merginggrids[k]==grid1) || (grids[k]==grid1))redundant=1;
				}
				if(!redundant){
					/*Ok, add grid1 to grids, and grid3 to merginggrids: */
					grids[counter]=grid1;
					merginggrids[counter]=grid3;
					counter++;
				}
				/*Check that grid2 is not already present in the merginggrids: */
				redundant=0;
				for (k=0;k<counter;k++){
					if ((merginggrids[k]==grid1) || (grids[k]==grid1))redundant=1;
				}
				if(!redundant){
					/*Ok, add grid2 to grids, and grid4 to merginggrids: */
					grids[counter]=grid2;
					merginggrids[counter]=grid4;
					counter++;
				}
			}
		}
	}

	/*Ok, we have counter pairs of grids (grids and merginggrids): start merging grids in the triangulation: */
	newnods=nods;
	for (i=0;i<counter;i++){
		grid1=grids[i];
		grid3=merginggrids[i];
		/*Merge grid3 into grid1: */ 
		x[grid3]=xmin; //flag  for later removal from x
		y[grid3]=ymin; //flag  for later removal from y
		newnods--;
		for (k=0;k<nel;k++){
			if (*(index+3*k+0)==grid3)*(index+3*k+0)=grid1;
			if (*(index+3*k+1)==grid3)*(index+3*k+1)=grid1;
			if (*(index+3*k+2)==grid3)*(index+3*k+2)=grid1;
		}
	}

	/*Reallocate x and y: */
	xreal=(double*)xrealloc(x,newnods*sizeof(double));
	yreal=(double*)xrealloc(y,newnods*sizeof(double));
	counter1=0;
	counter2=0;
	for (i=0;i<nods;i++){
		if (x[i]!=xmin){
			xreal[counter1]=x[i];
			counter1++;
		}
		if (y[i]!=ymin){
			yreal[counter2]=y[i];
			counter2++;
		}
	}
	xfree((void**)&x); x=xreal;
	xfree((void**)&y); y=yreal;

	/*Assign output pointers:*/
	*pindex=index;
	*px=x;
	*py=y;
	*pnods=newnods;
	*psegments=segments;
	*pnumsegs=numsegs;

	return noerr;
}

/******************************************************************************************************************************
                                   IsRiftPresent
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "IsRiftPresent"
	
int IsRiftPresent(int* priftflag,int* pnumrifts, double* segmentmarkerlist,int nsegs){

	int i;
	int noerr=1;
	
	/*output: */
	int riftflag=0;
	int numrifts=0;

	double maxmark=1; //default marker for regular segments

	/*Any marker >=2 indicates a certain rift: */
	numrifts=0;
	for (i=0;i<nsegs;i++){
		if (segmentmarkerlist[i]>maxmark){
			numrifts++;
			maxmark=segmentmarkerlist[i];
		}
	}
	if (numrifts)riftflag=1;

	/*Assign output pointers:*/
	*priftflag=riftflag;
	*pnumrifts=numrifts;

	return noerr;
}

/******************************************************************************************************************************
                                   OrderRifts
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "OrderRifts"
		
int OrderRifts(double** priftstips, double** riftssegments,double** riftspairs,int numrifts,int* riftsnumsegments,double* x,double* y){
	
	int noerr=1;
	int i,j,k,counter;

	/*intermediary: */
	double* riftsegments=NULL;
	double* riftpairs=NULL;
	int numsegs;

	/*ordering and copy: */
	int*    order=NULL;
	double* riftsegments_copy=NULL;
	double* riftpairs_copy=NULL;

	/*grid and element manipulation: */
	int     grid1,grid2,grid3,grid4,temp_grid,tip1,tip2,grid;
	double  el1,el2;
	int     already_ordered=0;

	/*output: */
	double* riftstips=NULL;


	/*Allocate byproduct of this routine, riftstips: */
	riftstips=(double*)xmalloc(numrifts*2*sizeof(double));

	/*Go through all rifts: */
	for (i=0;i<numrifts;i++){
		riftsegments=riftssegments[i];
		riftpairs=riftspairs[i];
		numsegs=riftsnumsegments[i];
	
			
		/*Allocate copy of riftsegments and riftpairs, 
		 *as well as ordering vector: */
		riftsegments_copy=(double*)xmalloc(numsegs*3*sizeof(double));
		riftpairs_copy=(double*)xmalloc(numsegs*2*sizeof(double));
		order=(int*)xmalloc(numsegs*sizeof(int));

		/*First find the tips, using the pairs. If a pair of elements has one grid in common, this grid is a rift tip: */
		tip1=-1;
		tip2=-1;

		for (j=0;j<numsegs;j++){
			el1=*(riftpairs+2*j+0);
			el2=*(riftpairs+2*j+1);
			grid1=(int)*(riftsegments+3*j+0);
			grid2=(int)*(riftsegments+3*j+1);
			/*Summary, el1 and el2 are facing one another across the rift. grid1 and grid2 belong to el1 and 
			 *are located on the rift. Find grid3 and grid4, grids belonging to el2 and located on the rift: */
			for (k=0;k<numsegs;k++){
				if (*(riftsegments+3*k+2)==el2){
					grid3=(int)*(riftsegments+3*k+0);
					grid4=(int)*(riftsegments+3*k+1);
					break;
				}
			}
			/* Make sure grid3 faces grid1 and grid4 faces grid2: */
			if ((x[grid1]==x[grid4]) && (y[grid1]==y[grid4])){
				/*Swap grid3 and grid4:*/
				temp_grid=grid3;
				grid3=grid4;
				grid4=temp_grid;
			}

			/*Figure out if a tip is on this element: */
			if (grid3==grid1){
				/*grid1 is a tip*/
				if (tip1==-1) {
					tip1=grid1;
					continue;
				}
				if ((tip2==-1) && (grid1!=tip1)){
					tip2=grid1;
					break;
				}
			}
		
			if (grid4==grid2){
				/*grid2 is a tip*/
				if (tip1==-1){
					tip1=grid2;
					continue;
				}
				if ((tip2==-1) && (grid2!=tip1)){
					tip2=grid2;
					break;
				}
			}
		}

		#ifdef _DEBUG_
			printf("Tips for rift#%i   (%i-%i)\n",i,tip1,tip2);
		#endif
		/*Record tips in riftstips: */
		*(riftstips+2*i+0)=(double)tip1;
		*(riftstips+2*i+1)=(double)tip2;


		/*We have the two tips for this rift.  Go from tip1 to tip2, and figure out the order in which segments are sequential. 
		 *Because two elements are connected to tip1, we chose one first, which defines the direction we are rotating along the rift. */
		grid=tip1;
		for (counter=0;counter<numsegs;counter++){
			for (j=0;j<numsegs;j++){
				grid1=(int)*(riftsegments+3*j+0);
				grid2=(int)*(riftsegments+3*j+1);
				
				if ((grid1==grid) || (grid2==grid)){
					/*Ok, this segment is connected to grid, plug its index into order, unless we already plugged it before: */
					already_ordered=0;
					for (k=0;k<counter;k++){
						if(order[k]==j){
							already_ordered=1;
							break;
						}
					}
					if (!already_ordered){
						order[counter]=j;
						if(grid1==grid){
							grid=grid2;
						}
						else if(grid2==grid){
							grid=grid1;
						}
						break;
					}
				}
			}
		}

		#ifdef _DEBUG_
		for (j=0;j<numsegs;j++){
			printf("%i\n",order[j]);
		}
		#endif

		/*Using the order vector, and the riftsegments_copy and riftspairs_copy, reorder the segments and the pairs: */
		for (j=0;j<numsegs;j++){
			*(riftsegments_copy+3*j+0)=*(riftsegments+3*order[j]+0);
			*(riftsegments_copy+3*j+1)=*(riftsegments+3*order[j]+1);
			*(riftsegments_copy+3*j+2)=*(riftsegments+3*order[j]+2);
			*(riftpairs_copy+2*j+0)=*(riftpairs+2*order[j]+0);
			*(riftpairs_copy+2*j+1)=*(riftpairs+2*order[j]+1);
		}
		
		for (j=0;j<numsegs;j++){
			*(riftsegments+3*j+0)=*(riftsegments_copy+3*j+0);
			*(riftsegments+3*j+1)=*(riftsegments_copy+3*j+1);
			*(riftsegments+3*j+2)=*(riftsegments_copy+3*j+2);
			*(riftpairs+2*j+0)=*(riftpairs_copy+2*j+0);
			*(riftpairs+2*j+1)=*(riftpairs_copy+2*j+1);
		}

		xfree((void**)&order);
		xfree((void**)&riftsegments_copy);
		xfree((void**)&riftpairs_copy);

	}

	/*Assign output pointer:*/
	*priftstips=riftstips;
	return noerr;
}

/******************************************************************************************************************************
                                   PenaltyPairs
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "PenaltyPairs"
			
int PenaltyPairs(double*** priftspenaltypairs,int** priftsnumpenaltypairs,int numrifts,double** riftssegments,
		int* riftsnumsegs,double** riftspairs,double* riftstips,double* x,double* y){
		

	int noerr=1;
	int i,j,k,k0;

	double el1,el2,grid1,grid2,grid3,grid4;
	double tip1,tip2,temp_grid;

	/*output: */
	double** riftspenaltypairs=NULL;
	double*  riftpenaltypairs=NULL;
	int*     riftsnumpenaltypairs=NULL;

	/*intermediary: */
	int numsegs;
	double* riftsegments=NULL;
	double* riftpairs=NULL;
	int counter;
	double normal[2];
	double length;
	int    k1,k2;

	/*Allocate: */
	riftspenaltypairs=(double**)xmalloc(numrifts*sizeof(double*));
	riftsnumpenaltypairs=(int*)xmalloc(numrifts*sizeof(int));

	for(i=0;i<numrifts;i++){
		numsegs=riftsnumsegs[i];
		riftsegments=riftssegments[i];
		riftpairs=riftspairs[i];
		tip1=*(riftstips+2*i+0);
		tip2=*(riftstips+2*i+1);

		/*allocate riftpenaltypairs, and riftnumpenaltypairs: */
		if((numsegs/2-1)!=0)riftpenaltypairs=(double*)xcalloc((numsegs/2-1)*7,sizeof(double));
		
		/*Go through only one flank of the rifts, not counting the tips: */
		counter=0;
		for(j=0;j<(numsegs/2);j++){
			el1=*(riftpairs+2*j+0);
			el2=*(riftpairs+2*j+1);
			grid1=*(riftsegments+3*j+0);
			grid2=*(riftsegments+3*j+1);
			/*Find segment index to recover grid3 and grid4, facing grid1 and grid2: */
			k0=-1;
			for(k=0;k<numsegs;k++){
				if(*(riftsegments+3*k+2)==el2){
					k0=k;
					break;
				}
			}
			grid3=*(riftsegments+3*k0+0);
			grid4=*(riftsegments+3*k0+1);

			/* Make sure grid3 faces grid1 and grid4 faces grid2: */
			if ((x[(int)grid1-1]==x[(int)grid4-1]) && (y[(int)grid1-1]==y[(int)grid4-1])){
				/*Swap grid3 and grid4:*/
				temp_grid=grid3;
				grid3=grid4;
				grid4=temp_grid;
			}	
			/*Ok, we have grid1 facing grid3, and grid2 facing grid4. Compute the normal to 
			 *this segment, and its length: */
			normal[0]=cos(atan2(x[(int)grid1-1]-x[(int)grid2-1],y[(int)grid2-1]-y[(int)grid1-1]));
			normal[1]=sin(atan2(x[(int)grid1-1]-x[(int)grid2-1],y[(int)grid2-1]-y[(int)grid1-1]));
			length=sqrt(pow(x[(int)grid2-1]-x[(int)grid1-1],(double)2)+pow(y[(int)grid2-1]-y[(int)grid1-1],(double)2));

			/*Be careful here, we want penalty loads on each grid, not on each segment. This means we cannot plug grid1,
			 * grid2, grid3 and grid4 directly into riftpenaltypairs. We need to include grid1, grid2, grid3 and grid4, 
			 * only once. We'll add the normals and the lengths : */

			if(grid1!=grid3){ //exclude tips from loads
				k1=-1;
				for(k=0;k<counter;k++){
					if( (*(riftpenaltypairs+k*7+0))==grid1){
						k1=k; 
						break;
					}
				}
				if(k1==-1){
					*(riftpenaltypairs+counter*7+0)=grid1;
					*(riftpenaltypairs+counter*7+1)=grid3;
					*(riftpenaltypairs+counter*7+2)=el1;
					*(riftpenaltypairs+counter*7+3)=el2;
					*(riftpenaltypairs+counter*7+4)=normal[0];
					*(riftpenaltypairs+counter*7+5)=normal[1];
					*(riftpenaltypairs+counter*7+6)=length/2;
					counter++;
				}
				else{
					*(riftpenaltypairs+k1*7+4)+=normal[0];
					*(riftpenaltypairs+k1*7+5)+=normal[1];
					*(riftpenaltypairs+k1*7+6)+=length/2;
				}
			}
			if(grid2!=grid4){
				k2=-1;
				for(k=0;k<counter;k++){
					if( (*(riftpenaltypairs+k*7+0))==grid2){
						k2=k;
						break;
					}
				}
				if(k2==-1){
					*(riftpenaltypairs+counter*7+0)=grid2;
					*(riftpenaltypairs+counter*7+1)=grid4;
					*(riftpenaltypairs+counter*7+2)=el1;
					*(riftpenaltypairs+counter*7+3)=el2;
					*(riftpenaltypairs+counter*7+4)=normal[0];
					*(riftpenaltypairs+counter*7+5)=normal[1];
					*(riftpenaltypairs+counter*7+6)=length/2;
					counter++;
				}
				else{
					*(riftpenaltypairs+k2*7+4)+=normal[0];
					*(riftpenaltypairs+k2*7+5)+=normal[1];
					*(riftpenaltypairs+k2*7+6)+=length/2;
				}
			}
		}
		/*Renormalize normals: */
		for(j=0;j<counter;j++){
			double magnitude=sqrt(pow( *(riftpenaltypairs+j*7+4),2) + pow( *(riftpenaltypairs+j*7+5),2) );
			*(riftpenaltypairs+j*7+4)=*(riftpenaltypairs+j*7+4)/magnitude;
			*(riftpenaltypairs+j*7+5)=*(riftpenaltypairs+j*7+5)/magnitude;
		}
		
		riftspenaltypairs[i]=riftpenaltypairs;
		riftsnumpenaltypairs[i]=(numsegs/2-1);
	}
			


	/*Assign output pointers: */
	*priftspenaltypairs=riftspenaltypairs;
	*priftsnumpenaltypairs=riftsnumpenaltypairs;

	return noerr;
}

/******************************************************************************************************************************
                                   RemoveCorners
******************************************************************************************************************************/

#undef __FUNCT__ 
#define __FUNCT__ "RemoveCorners"
		

int RemoveCornersFromRifts(double** pindex,int* pnel,double** px,double** py,int* pnods, double* segments,double* segmentmarkers,int num_seg){

	int noerr=1;
	int i,j,k;
	double grid1,grid2,grid3;
	int el;

	/*input: */
	double* index=NULL;
	int     nel;
	double* x=NULL;
	double* y=NULL;
	int     nods;
	double  pair[2];
	int     pair_count=0;
	int     triple=0;

	/*Recover input: */
	index=*pindex;
	nel=*pnel;
	x=*px;
	y=*py;
	nods=*pnods;


	for (i=0;i<num_seg;i++){
		grid1=*(segments+3*i+0);
		grid2=*(segments+3*i+1);
		/*Find all elements connected to [grid1 grid2]: */
		pair_count=0;
		for (j=0;j<nel;j++){
			if (*(index+3*j+0)==grid1){
				if ((*(index+3*j+1)==grid2) || (*(index+3*j+2)==grid2)){
					pair[pair_count]=j;
					pair_count++;
				}
			}
			if (*(index+3*j+1)==grid1){
				if ((*(index+3*j+0)==grid2) || (*(index+3*j+2)==grid2)){
					pair[pair_count]=j;
					pair_count++;
				}
			}
			if (*(index+3*j+2)==grid1){
				if ((*(index+3*j+0)==grid2) || (*(index+3*j+1)==grid2)){
					pair[pair_count]=j;
					pair_count++;
				}
			}
		}
		/*Ok, we have pair_count elements connected to this segment. For each of these elements, 
		 *figure out if the third grid also belongs to a segment: */
		if ((pair_count==0) || (pair_count==1)){ //we only select the rift segments, which belong to  2 elements
			continue;
		}
		else{
			for (j=0;j<pair_count;j++){
				el=(int)pair[j];
				triple=0;
				/*First find grid3: */
				if (*(index+3*el+0)==grid1){
					if (*(index+3*el+1)==grid2)grid3=*(index+3*el+2);
					else grid3=*(index+3*el+1);
				}
				if (*(index+3*el+1)==grid1){
					if (*(index+3*el+0)==grid2)grid3=*(index+3*el+2);
					else grid3=*(index+3*el+0);
				}
				if (*(index+3*el+2)==grid1){
					if (*(index+3*el+0)==grid2)grid3=*(index+3*el+1);
					else grid3=*(index+3*el+0);
				}
				/*Ok, we have grid3. Does grid3 belong to a segment? : */
				for (k=0;k<num_seg;k++){
					if ((grid3==*(segments+3*k+0)) || (grid3==*(segments+3*k+1))){
						triple=1;
						break;
					}
				}
				if(triple==1){
					/*el is a corner element: we need to split it in 3 triangles: */
					x=(double*)xrealloc(x,(nods+1)*sizeof(double));
					y=(double*)xrealloc(y,(nods+1)*sizeof(double));
					x[nods]=(x[(int)grid1-1]+x[(int)grid2-1]+x[(int)grid3-1])/3;
					y[nods]=(y[(int)grid1-1]+y[(int)grid2-1]+y[(int)grid3-1])/3;

					index=(double*)xrealloc(index,(nel+2)*3*sizeof(double));
					/*First, reassign element el: */
					*(index+3*el+0)=grid1;
					*(index+3*el+1)=grid2;
					*(index+3*el+2)=nods+1;
					/*Other two elements: */
					*(index+3*nel+0)=grid2;
					*(index+3*nel+1)=grid3;
					*(index+3*nel+2)=nods+1;

					*(index+3*(nel+1)+0)=grid3;
					*(index+3*(nel+1)+1)=grid1;
					*(index+3*(nel+1)+2)=nods+1;
					/*we need  to change the segment elements corresponding to el: */
					for (k=0;k<num_seg;k++){
						if (*(segments+3*k+2)==(double)(el+1)){
							if ( ((*(segments+3*k+0)==grid1) && (*(segments+3*k+1)==grid2)) || ((*(segments+3*k+0)==grid2) && (*(segments+3*k+1)==grid1))) *(segments+3*k+2)=(double)(el+1);
							if ( ((*(segments+3*k+0)==grid2) && (*(segments+3*k+1)==grid3)) || ((*(segments+3*k+0)==grid3) && (*(segments+3*k+1)==grid2))) *(segments+3*k+2)=(double)(nel+1);
							if ( ((*(segments+3*k+0)==grid3) && (*(segments+3*k+1)==grid1)) || ((*(segments+3*k+0)==grid1) && (*(segments+3*k+1)==grid3))) *(segments+3*k+2)=(double)(nel+2);
						}
					}

					nods=nods+1;
					nel=nel+2;
					i=0;
					break;
				}
			} //for (j=0;j<pair_count;j++)
		}
	}// for (i=0;i<num_seg;i++)

	/*Assign output pointers: */
	*pindex=index;
	*pnel=nel;
	*px=x;
	*py=y;
	*pnods=nods;

	return noerr;
}
