/*!\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[10];
	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;

	/* state: */
	double* state=NULL;

	int     num_seg;

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

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

	/*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{
		_error_("first argument should be the element list");
	}

	/*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{
		_error_("second argument should be the x corrdinate list");
	}

	/*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{
		_error_("third argument should be the y corrdinate list");
	}	

	/*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{
		_error_("fourth argument should be the segments list");
	}

	/*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{
		_error_("fourth argument should be the segmentmarkers list");
	}

	/*First, do some fixing on the existing mesh: we do not want any element belonging entirely to the segment list (ie: 
	 *all the nodes 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: */
	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,nods,nel);

		/*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,nods,nel);

		/*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 : */
	WriteData(&plhs[0],index_in,nel,3);
	//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, "transpose");

	WriteData(&plhs[1],x_in,nods,1);
	//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, "transpose");

	WriteData(&plhs[2],y_in,nods,1);
	//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, "transpose");

	WriteData(&plhs[3],segments_in,num_seg,3);
//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, "transpose");

	WriteData(&plhs[4],segmentmarkers_in,num_seg,1);
	//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, "transpose");

	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";
		fnames[7] = "fraction";
		fnames[8] = "fractionincrement";
		fnames[9] = "state";

		dimensions[0]=out_numrifts;

		pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
		
		for (i=0;i<out_numrifts;i++){

			/*Segments: */
			WriteData(&pmxa_array3,out_riftssegments[i],out_riftsnumsegments[i],3);
			//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, "transpose");
			mxSetField(pmxa_array,i,"segments",pmxa_array3);
			mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)out_riftsnumsegments[i]));

			/*Element pairs: */
			WriteData(&pmxa_array3,out_riftspairs[i],out_riftsnumpairs[i],2);
			//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, "transpose");
			mxSetField(pmxa_array,i,"pairs",pmxa_array3);

			/*Tips: */
			pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
			mxSetM(pmxa_array2,1);
			pair=(double*)mxMalloc(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: */
			WriteData(&pmxa_array3,out_riftspenaltypairs[i],out_riftsnumpenaltypairs[i],7);
			//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, "transpose");
			mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);

			/*Friction fraction, fractionincrement  and fill: */
			mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
			mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum)); //default is ice
			mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
			mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1)); 

			/*State: */
			state=(double*)mxMalloc(out_riftsnumpenaltypairs[i]*sizeof(double));
			for(j=0;j<out_riftsnumpenaltypairs[i];j++)state[j]=FreeEnum;
			pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
			mxSetM(pmxa_array2,1);
			mxSetN(pmxa_array2,out_riftsnumpenaltypairs[i]);
			mxSetPr(pmxa_array2,state);
			mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
			
			mxSetField(pmxa_array,i,"state",pmxa_array3);
		}
	}
	else{
		/*output NaN :*/
		pNaN=(double*)mxMalloc(sizeof(double));
		*pNaN=NAN;
		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)
{
	_printLine_("");
	_printLine_("   usage: [index2,x2,y2,segments2,segmentmarkers2,rifts2]=TriMeshProcessrifts(index1,x1,y1,segments1,segmentmarkers1) ");
	_printLine_("      where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation.");
	_printLine_("      index2,x2,y2,segments2,segmentmarkers2,rifts2 is the resulting triangulation where rifts have been processed.");
}