/*!\file:  TriMeshRefine
 * \brief refine a mesh output by TriMesh:
 */ 

#include "./TriMeshRefine.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: */

	double* index=NULL;
	double* x=NULL;
	double* y=NULL;
	double* segments=NULL;
	double* segmentmarkerlist=NULL;
	
	/* input: */
	double* area=NULL;
	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;
	char*   order=NULL;

	/* error handling: */
	int		noerr=1;	
  
	/* Triangle structures: */
	struct triangulateio in,out;
	char   options[256];

	/* verify correct usage: */
	if (nlhs==0 && nrhs==0) {
		/* special case: */
		TriMeshRefineUsage();
		return;
	}

	if (!(  (nlhs==5) || (nrhs==6) || (nrhs==7)  )){
		mexPrintf("   %s format error.\n", __FUNCT__);
		TriMeshRefineUsage();
		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)-1;
			}
		}
	}
	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)-1;
			}
		}
	}
	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(" ");
	}
	
	/*Fetch area: */
	if(mxIsDouble(prhs[5])){
		if (mxGetM(prhs[5])!=nel){
			printf("%s%s\n",__FUNCT__," error message: area vector should be of the same size as index\n");
			mexErrMsgTxt(" ");
		}
		area=mxGetPr(prhs[5]);
	}
	else{
		printf("%s%s\n",__FUNCT__," error message: fifth argument should be the area list!");
		mexErrMsgTxt(" ");
	}

	/*Optionaly, recover desired order: */
	if(nrhs==7){
		if (!mxIsChar(prhs[6])){
			mexPrintf("%s%s\n",__FUNCT__," error message; sixth argument should be a string ('yes' or 'no')!");
			mexErrMsgTxt(" ");
		}
		order = (char *) xmalloc((mxGetN(prhs[6])+1)*sizeof(char));
		mxGetString(prhs[6],order,mxGetN(prhs[6])+1);
	}


	#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");
	}*/
	/*printf("area: \n");
	for (i=0;i<nel;i++){
		printf("%16.16lf \n",area[i]);
	}*/
	//if(nrhs==6)printf("Order: %s\n",order);
	#endif

	/*Create initial triangulation to call triangulate():*/
	in.numberoftriangles=nel;
	in.numberoftriangleattributes=1;
	in.numberofcorners=3;
	
	in.trianglelist = (int *) xmalloc(3*in.numberoftriangles * sizeof(int));
	for(i=0;i<in.numberoftriangles;i++){
		for(j=0;j<3;j++){
			in.trianglelist[3*i+j]=(int)index_in[3*i+j];
		}
	}
	in.triangleattributelist = (REAL *) xmalloc(in.numberoftriangles * in.numberoftriangleattributes * sizeof(REAL));
	for(i=0;i<in.numberoftriangles;i++){
		in.triangleattributelist[i]=0.0;
	}
	in.trianglearealist = (REAL *) xmalloc(in.numberoftriangles * sizeof(REAL));
	for(i=0;i<in.numberoftriangles;i++){
		in.trianglearealist[i]=area[i];
	}

	in.numberofpoints=nods;
	in.numberofpointattributes=1;
	in.pointlist = (REAL *) xmalloc(in.numberofpoints * 2 * sizeof(REAL));
	for (i=0;i<nods;i++){
		in.pointlist[2*i+0]=x_in[i];
		in.pointlist[2*i+1]=y_in[i];
	}
	in.pointattributelist = (REAL *) xmalloc(in.numberofpoints * in.numberofpointattributes * sizeof(REAL));
	for (i=0;i<nods;i++){
		in.pointattributelist[i] = 0.0;
	}	

	in.numberofsegments = num_seg;
	in.segmentlist = (int *) xmalloc(in.numberofsegments * 2 * sizeof(REAL));
	in.segmentmarkerlist = (int *) mxCalloc(in.numberofsegments,sizeof(int));
	for (i=0;i<num_seg;i++){
		in.segmentlist[2*i+0]=(int)segments_in[3*i+0];
		in.segmentlist[2*i+1]=(int)segments_in[3*i+1];
		in.segmentmarkerlist[i]=(int)segmentmarkers_in[i];
	}
	
	in.numberofholes = 0;

	/* Make necessary initializations so that Triangle can return a */
	/*   triangulation in `out'.  */

	out.pointlist = (REAL *) NULL;            
	out.pointattributelist = (REAL *) NULL;
	out.pointmarkerlist = (int *) NULL; 
	out.trianglelist = (int *) NULL;          
	out.triangleattributelist = (REAL *) NULL;
	out.neighborlist = (int *) NULL;         
	out.segmentlist = (int *) NULL;
	out.segmentmarkerlist = (int *) NULL;
	out.edgelist = (int *) NULL;             
	out.edgemarkerlist = (int *) NULL;   

	/* Triangulate the points:.  Switches are chosen to read and write a  */
	/*   PSLG (p), preserve the convex hull (c), number everything from  */
	/*   zero (z), assign a regional attribute to each element (A), and  */
	/*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
	/*   neighbor list (n).                                              */

	sprintf(options,"%s%lf","QzDq30iarp",area); //replace V by Q to quiet down the logging


	triangulate(options, &in, &out, NULL);
	
	/*Allocate index, x and y: */
	index=(double*)xmalloc(3*out.numberoftriangles*sizeof(double));
	x=(double*)xmalloc(out.numberofpoints*sizeof(double));
	y=(double*)xmalloc(out.numberofpoints*sizeof(double));
	segments=(double*)xmalloc(3*out.numberofsegments*sizeof(double));
	segmentmarkerlist=(double*)xmalloc(out.numberofsegments*sizeof(double));

	for (i = 0; i < out.numberoftriangles; i++) {
		for (j = 0; j < out.numberofcorners; j++) {
			*(index+3*i+j)=(double)out.trianglelist[i * out.numberofcorners + j]+1;
		}
	}
	for (i = 0; i < out.numberofpoints; i++) {
		x[i]=out.pointlist[i * 2 + 0];
		y[i]=out.pointlist[i * 2 + 1];
	}
	
	for (i = 0; i < out.numberofsegments; i++) {
		segments[3*i+0]=(double)out.segmentlist[i*2+0]+1;
		segments[3*i+1]=(double)out.segmentlist[i*2+1]+1;
		segmentmarkerlist[i]=(double)out.segmentmarkerlist[i];
	}

	/*Associate elements with segments: */
	noerr=AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);

	/*Order segments so that their normals point outside the domain: */
	if(!strcmp(order,"yes")){
		noerr=OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
	}
	
	/*Output : */
	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,3);
	mxSetN(pmxa_array,out.numberoftriangles);
	mxSetPr(pmxa_array,index);
	mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,out.numberofpoints);
	mxSetPr(pmxa_array,x);
	mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,out.numberofpoints);
	mxSetPr(pmxa_array,y);
	mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,3);
	mxSetN(pmxa_array,out.numberofsegments);
	mxSetPr(pmxa_array,segments);
	mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "'");

	pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetM(pmxa_array,1);
	mxSetN(pmxa_array,out.numberofsegments);
	mxSetPr(pmxa_array,segmentmarkerlist);
	mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "'");

	return;
}


void TriMeshRefineUsage(void)
{
	printf("\n");
	printf("   usage: [index2,x2,y2,segments2,segmentmarkers2]=TriMeshRefine(index1,x1,y1,segments1,segmentmarkers1,area,order) \n");
	printf("      where: (index1,x1,y1,segments1,segmentmarkers1) is an initial triangulation (segment1 does not need to be ordered).\n");
	printf("      area is a vector of element areas determining which elements will be refined.\n");
	printf("      order is an optional argument that determines whether segments are output in the \n");
	printf("      order they are made by Triangle (ie none), or ordered counter clockwise around the domain outline.\n");
	printf("      (index2,x2,y2,segments2,segmentmarkers2) is the resulting refined triangulation.\n");
	printf("   note: order is an optional arguments\n");
}
