/*!\file:  TriMeshProcessRifts.cpp
 * \brief split a mesh where a rift (or fault) is present
 */ 

#include "./TriMeshProcessRifts.h"

void mexFunction(	int nlhs, mxArray* plhs[],
					int nrhs, const mxArray* prhs[] ) {


	/*Matlab arrays: */
	mxArray* pmxa_array=NULL;
	mxArray* pmxa_array2=NULL;
	mxArray* pmxa_array3=NULL;
	int i,j,k,counter;
	
	/* returned quantities: */
	int      out_numrifts;
	int*     out_riftsnumsegments=NULL;
	double** out_riftssegments=NULL; 
	int*     out_riftsnumpairs=NULL;
	double** out_riftspairs=NULL;
	double*  out_riftstips=NULL;
	double** out_riftspenaltypairs=NULL;
	int*     out_riftsnumpenaltypairs=NULL;

	/*empty rifts structure: */
	double* pNaN=NULL;
	const	char*	fnames[7];
	mwSize     ndim=2;
	mwSize		dimensions[2] = {1,1};
	double* pair=NULL;

	
	/* input: */
	double* tindex_in=NULL;
	double* index_in=NULL;
	int     nel;
	double* x_inm=NULL; //matlab vector
	double* x_in=NULL; //copy of matlab vector
	int     nods;
	double* y_inm=NULL;//matlab vector
	double* y_in=NULL;//copy of matlab vector
	double* tsegments_in=NULL;
	double* segments_in=NULL;
	double* tsegmentmarkers_in=NULL;
	double* segmentmarkers_in=NULL;

	int     num_seg;

	/*rifts: */
	int     riftflag;
	int     numrifts;

	/* error handling: */
	int		noerr=1;	
  
	/* verify correct usage: */
	if (nlhs==0 && nrhs==0) {
		/* special case: */
		TriMeshProcessRiftsUsage();
		return;
	}
	
	if (!(  (nlhs==6) || (nrhs==5))){
		mexPrintf("   %s format error.\n", __FUNCT__);
		TriMeshProcessRiftsUsage();
		printf("   ");
		mexErrMsgTxt(" ");
	}

	/*Fetch index_in: */
	if(mxIsDouble(prhs[0])){
		nel=mxGetM(prhs[0]);
		tindex_in=mxGetPr(prhs[0]);
		index_in=(double*)xmalloc(nel*3*sizeof(double));
		for (i=0;i<nel;i++){
			for (j=0;j<3;j++){
				*(index_in+3*i+j)=*(tindex_in+nel*j+i);
			}
		}
	}
	else{
		printf("%s%s\n",__FUNCT__," error message: first argument should be the element list!");
		mexErrMsgTxt(" ");
	}

	/*Fetch x_in: */
	if(mxIsDouble(prhs[1])){
		nods=mxGetM(prhs[1]);
		x_inm=mxGetPr(prhs[1]);
		x_in=(double*)xmalloc(nods*sizeof(double));
		for (i=0;i<nods;i++){
			x_in[i]=x_inm[i];
		}
	}
	else{
		printf("%s%s\n",__FUNCT__," error message: second argument should be the x corrdinate list!");
		mexErrMsgTxt(" ");
	}

	/*Fetch y_in: */
	if(mxIsDouble(prhs[2])){
		y_inm=mxGetPr(prhs[2]);
		y_in=(double*)xmalloc(nods*sizeof(double));
		for (i=0;i<nods;i++){
			y_in[i]=y_inm[i];
		}
	}
	else{
		printf("%s%s\n",__FUNCT__," error message: third argument should be the y corrdinate list!");
		mexErrMsgTxt(" ");
	}	

	/*Fetch segments_in: */
	if(mxIsDouble(prhs[3])){
		num_seg=mxGetM(prhs[3]);
		tsegments_in=mxGetPr(prhs[3]);
		segments_in=(double*)xmalloc(num_seg*3*sizeof(double));
		for (i=0;i<num_seg;i++){
			for (j=0;j<3;j++){
				*(segments_in+3*i+j)=*(tsegments_in+num_seg*j+i);
			}
		}
	}
	else{
		printf("%s%s\n",__FUNCT__," error message: fourth argument should be the segments list!");
		mexErrMsgTxt(" ");
	}

	/*Fetch segment markers: */
	if(mxIsDouble(prhs[4])){
		tsegmentmarkers_in=mxGetPr(prhs[4]);
		segmentmarkers_in=(double*)xmalloc(num_seg*sizeof(double));
		for (i=0;i<num_seg;i++){
			segmentmarkers_in[i]=tsegmentmarkers_in[i];
		}
	}
	else{
		printf("%s%s\n",__FUNCT__," error message: fourth argument should be the segmentmarkers list!");
		mexErrMsgTxt(" ");
	}

	#ifdef _DEBUG_
	printf("Index: \n");
	for (i=0;i<nel;i++){
		for(j=0;j<3;j++){
			printf("%lf ",*(index_in+3*i+j));
		}
		printf("\n");
	}
	printf("x,y: \n");
	for (i=0;i<nods;i++){
		printf("%16.16lf %16.16lf\n",x_in[i],y_in[i]);
	}
	printf("segments:\n");
	for (i=0;i<num_seg;i++){
		for(j=0;j<3;j++){
			printf("%lf ",*(segments_in+3*i+j));
		}
		printf("%lf ",segmentmarkers_in[i]);
		printf("\n");
	}
	#endif

	/*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie: 
	 *all the grids of this element belong to the segments (tends to happen when there are corners: */
	RemoveCornersFromRifts(&index_in,&nel,&x_in,&y_in,&nods,segments_in,segmentmarkers_in,num_seg);

	/*Figure out if we have rifts, and how many: */
	IsRiftPresent(&riftflag,&numrifts,segmentmarkers_in,num_seg);
	
	if(riftflag){	
		SplitMeshForRifts(&nel,&index_in,&nods,&x_in,&y_in,&num_seg,&segments_in,&segmentmarkers_in);
	}

	/*Order segments so that their normals point outside the domain: */
	noerr=OrderSegments(&segments_in,num_seg, index_in,nel);

	
	if(riftflag){
		
		/*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the 
		 *segmentmarkerlist:*/
		SplitRiftSegments(&segments_in,&segmentmarkers_in,&num_seg,&out_numrifts,&out_riftsnumsegments,&out_riftssegments,numrifts);

		/*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
		PairRiftElements(&out_riftsnumpairs,&out_riftspairs,out_numrifts,out_riftsnumsegments,out_riftssegments,x_in,y_in);
		
		/*Order rifts so that they start from one tip, go to the other tip, and back: */
		OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in);

		/*Create penalty pairs, used by Imp: */
		PenaltyPairs(&out_riftspenaltypairs,&out_riftsnumpenaltypairs,numrifts,out_riftssegments,out_riftsnumsegments,out_riftspairs,out_riftstips,x_in,y_in);
	}


	/*Output : */
	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,3);
	mxSetN(pmxa_array,nel);
	mxSetPr(pmxa_array,index_in);
	mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,nods);
	mxSetPr(pmxa_array,x_in);
	mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,nods);
	mxSetPr(pmxa_array,y_in);
	mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "'");

	
	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,3);
	mxSetN(pmxa_array,num_seg);
	mxSetPr(pmxa_array,segments_in);
	mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "'");
	

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,num_seg);
	mxSetPr(pmxa_array,segmentmarkers_in);
	mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "'");

	if(riftflag){
		/*Create a structure rifts where if i is a rift number, we have the following fields rifts(i).segments and rifts(i).numsegs: */

		fnames[0] = "numsegs";
		fnames[1] = "segments";
		fnames[2] = "pairs";
		fnames[3] = "tips";
		fnames[4] = "penaltypairs";
		fnames[5] = "fill";
		fnames[6] = "friction";

		dimensions[0]=out_numrifts;

		pmxa_array=mxCreateStructArray( ndim,dimensions,7,fnames);
		
		for (i=0;i<out_numrifts;i++){
			/*Segments: */
			pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
			mxSetM(pmxa_array2,3);
			mxSetN(pmxa_array2,out_riftsnumsegments[i]);
			mxSetPr(pmxa_array2,out_riftssegments[i]);
			mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "'");
			
			mxSetField(pmxa_array,i,"segments",pmxa_array3);
			mxSetField(pmxa_array,i,"numsegs",mxCreateScalarDouble((double)out_riftsnumsegments[i]));

			/*Element pairs: */
			pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
			mxSetM(pmxa_array2,2);
			mxSetN(pmxa_array2,out_riftsnumpairs[i]);
			mxSetPr(pmxa_array2,out_riftspairs[i]);
			mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "'");
			
			mxSetField(pmxa_array,i,"pairs",pmxa_array3);

			/*Tips: */
			pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
			mxSetM(pmxa_array2,1);
			pair=(double*)xmalloc(2*sizeof(double));
			pair[0]=*(out_riftstips+2*i+0);
			pair[1]=*(out_riftstips+2*i+1);
			mxSetN(pmxa_array2,2);
			mxSetPr(pmxa_array2,pair);
			mxSetField(pmxa_array,i,"tips",pmxa_array2);

			/*Penalty pairs: */
			pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
			mxSetM(pmxa_array2,7);
			mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
			mxSetPr(pmxa_array2,out_riftspenaltypairs[i]);
			mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "'");
			
			mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);

			/*Friction and fill: set to 0 both */
			mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
			mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(1)); //default is water



		}
	}
	else{
		/*output NaN :*/
		pNaN=(double*)xmalloc(sizeof(double));
		*pNaN=0.0/0.0;
		pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
		mxSetM(pmxa_array,1);
		mxSetN(pmxa_array,1);
		mxSetPr(pmxa_array,pNaN);
		
	}
	plhs[5]=pmxa_array;

	return;
}


void TriMeshProcessRiftsUsage(void)
{
	printf("\n");
	printf("   usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(index1,x1,y1,segments1,segmentmarkers1) \n");
	printf("      where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.\n");
	printf("      index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.\n");
}
