Index: ../trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp
===================================================================
--- ../trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp	(revision 13587)
+++ ../trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp	(revision 13588)
@@ -10,48 +10,31 @@
 	_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.");
 }/*}}}*/
+void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg,
+			int  *pnumrifts, int  **priftsnumsegments, double ***priftssegments,
+			int  **priftsnumpairs, double ***priftspairs, double **priftstips,
+			double ***priftspenaltypairs, int  **priftsnumpenaltypairs);
 WRAPPER(TriMeshProcessRifts){
 
-	/*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;
+	int      numrifts;
+	int     *riftsnumsegments     = NULL;
+	double **riftssegments        = NULL;
+	int     *riftsnumpairs        = NULL;
+	double **riftspairs           = NULL;
+	double  *riftstips            = NULL;
+	double **riftspenaltypairs    = NULL;
+	int     *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 *index_in  = NULL;
-	int     nel;
-	double *x_in  = NULL; //copy of matlab vector
-	int     nods;
-	double *y_in               = NULL; //copy of matlab vector
-	double *segments_in        = NULL;
-	double *segmentmarkers_in  = NULL;
-
-	/* state: */
-	double *state = NULL;
+	int     nel,nods;
+	double *index          = NULL;
+	double *x              = NULL;
+	double *y              = NULL;
+	double *segments       = NULL;
+	double *segmentmarkers = NULL;
 	int     num_seg;
 
-	/*rifts: */
-	int     riftflag;
-	int     numrifts;
-
 	/*Boot module*/
 	MODULEBOOT();
 
@@ -59,121 +42,168 @@
 	CHECKARGUMENTS(NLHS,NRHS,&TriMeshProcessRiftsUsage);
 
 	/*Fetch data */
-	FetchData(&index_in,&nel,NULL,INDEXIN);
-	FetchData(&x_in,&nods,XIN);
-	FetchData(&y_in,NULL,YIN);
-	FetchData(&segments_in,&num_seg,NULL,SEGMENTSIN);
-	FetchData(&segmentmarkers_in,NULL,SEGMENTMARKERSIN);
+	FetchData(&index,&nel,NULL,INDEXIN);
+	FetchData(&x,&nods,XIN);
+	FetchData(&y,NULL,YIN);
+	FetchData(&segments,&num_seg,NULL,SEGMENTSIN);
+	FetchData(&segmentmarkers,NULL,SEGMENTMARKERSIN);
 
-	/*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);
+	/*call x layer*/
+	TriMeshProcessRiftsx(&index,&nel,&x,&y,&nods,&segments,&segmentmarkers,&num_seg,
+				&numrifts,&riftsnumsegments,&riftssegments,&riftsnumpairs,&riftspairs,&riftstips,&riftspenaltypairs,&riftsnumpenaltypairs);
 
-	/*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);
-	}
+	/*Output : */
+	WriteData(INDEXOUT,index,nel,3);
+	WriteData(XOUT,x,nods,1);
+	WriteData(YOUT,y,nods,1);
+	WriteData(SEGMENTSOUT,segments,num_seg,3);
+	WriteData(SEGMENTMARKERSOUT,segmentmarkers,num_seg,1);
 
-	/*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);
+	mxArray *pmxa_array  = NULL;
+	mxArray *pmxa_array2 = NULL;
+	mxArray *pmxa_array3 = NULL;
+	const char *fnames[10];
+	mwSize      ndim          = 2;
+	mwSize      dimensions[2] = {1,1};
+	double     *pair          = NULL;
 
-		/*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);
-	}
+	/*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]=numrifts;
+	pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
 
-	/*Output : */
-	WriteData(INDEXOUT,index_in,nel,3);
-	WriteData(XOUT,x_in,nods,1);
-	WriteData(YOUT,y_in,nods,1);
-	WriteData(SEGMENTSOUT,segments_in,num_seg,3);
-	WriteData(SEGMENTMARKERSOUT,segmentmarkers_in,num_seg,1);
+	for(int i=0;i<numrifts;i++){
 
-	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";
+		/*Segments: */
+		WriteData(&pmxa_array3,riftssegments[i],riftsnumsegments[i],3);
+		mxSetField(pmxa_array,i,"segments",pmxa_array3);
+		mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)riftsnumsegments[i]));
 
-		dimensions[0]=out_numrifts;
-		pmxa_array=mxCreateStructArray( ndim,dimensions,10,fnames);
-		
-		for (i=0;i<out_numrifts;i++){
+		/*Element pairs: */
+		WriteData(&pmxa_array3,riftspairs[i],riftsnumpairs[i],2);
+		mxSetField(pmxa_array,i,"pairs",pmxa_array3);
 
-			/*Segments: */
-			WriteData(&pmxa_array3,out_riftssegments[i],out_riftsnumsegments[i],3);
-			mxSetField(pmxa_array,i,"segments",pmxa_array3);
-			mxSetField(pmxa_array,i,"numsegs",mxCreateDoubleScalar((double)out_riftsnumsegments[i]));
+		/*Tips: */
+		pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
+		mxSetM(pmxa_array2,1);
+		pair=(double*)mxMalloc(2*sizeof(double));
+		pair[0]=*(riftstips+2*i+0);
+		pair[1]=*(riftstips+2*i+1);
+		mxSetN(pmxa_array2,2);
+		mxSetPr(pmxa_array2,pair);
+		mxSetField(pmxa_array,i,"tips",pmxa_array2);
 
-			/*Element pairs: */
-			WriteData(&pmxa_array3,out_riftspairs[i],out_riftsnumpairs[i],2);
-			mxSetField(pmxa_array,i,"pairs",pmxa_array3);
+		/*Penalty pairs: */
+		WriteData(&pmxa_array3,riftspenaltypairs[i],riftsnumpenaltypairs[i],7);
+		mxSetField(pmxa_array,i,"penaltypairs",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);
+		/*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)); 
 
-			/*Penalty pairs: */
-			WriteData(&pmxa_array3,out_riftspenaltypairs[i],out_riftsnumpenaltypairs[i],7);
-			mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
+		/*State: */
+		double *state=(double*)mxMalloc(riftsnumpenaltypairs[i]*sizeof(double));
+		for(int j=0;j<riftsnumpenaltypairs[i];j++) state[j]=FreeEnum;
+		pmxa_array2= mxCreateDoubleMatrix(0,0,mxREAL);
+		mxSetM(pmxa_array2,1);
+		mxSetN(pmxa_array2,riftsnumpenaltypairs[i]);
+		mxSetPr(pmxa_array2,state);
+		mexCallMATLAB( 1, &pmxa_array3, 1, &pmxa_array2, "transpose");
 
-			/*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)); 
+		mxSetField(pmxa_array,i,"state",pmxa_array3);
+	}
 
-			/*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;
 
 	/*end module: */
 	MODULEEND();
 }
+void TriMeshProcessRiftsx(double** pindex, int* pnel,double** px,double** py,int* pnods,double** psegments,double** psegmentmarkers,int *pnum_seg,
+			int  *pnumrifts, int  **priftsnumsegments, double ***priftssegments,
+			int  **priftsnumpairs, double ***priftspairs, double **priftstips,
+			double ***priftspenaltypairs, int  **priftsnumpenaltypairs){
+
+	/*Output*/
+	int      numrifts,numrifts0;
+	int     *riftsnumsegments     = NULL;
+	double **riftssegments        = NULL;
+	int     *riftsnumpairs        = NULL;
+	double **riftspairs           = NULL;
+	double  *riftstips            = NULL;
+	double **riftspenaltypairs    = NULL;
+	int     *riftsnumpenaltypairs = NULL;
+
+	/*Recover initial mesh*/
+	int     nel            = *pnel;
+	double *index          = *pindex;
+	double *x              = *px;
+	double *y              = *py;
+	int     nods           = *pnods;
+	double *segments       = *psegments;
+	double *segmentmarkers = *psegmentmarkers;
+	int     num_seg        = *pnum_seg;
+
+	/*Intermediary*/
+	int     riftflag;
+
+	/*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,&nel,&x,&y,&nods,segments,segmentmarkers,num_seg);
+
+	/*Figure out if we have rifts, and how many: */
+	IsRiftPresent(&riftflag,&numrifts0,segmentmarkers,num_seg);
+
+	if(!riftflag) _error_("No rift present in mesh");
+
+	/*Split mesh*/
+	SplitMeshForRifts(&nel,&index,&nods,&x,&y,&num_seg,&segments,&segmentmarkers);
+
+	/*Order segments so that their normals point outside the domain: */
+	OrderSegments(&segments,num_seg, index,nel);
+
+	/*We do not want to output segments mixed with rift segments: wring out the rifts from the segments, using the 
+	 *segmentmarkerlist:*/
+	SplitRiftSegments(&segments,&segmentmarkers,&num_seg,&numrifts,&riftsnumsegments,&riftssegments,numrifts0,nods,nel);
+
+	/*Using rift segments, associate rift faces in pairs, each pair face representing opposite flanks of the rifts facing one another directly: */
+	PairRiftElements(&riftsnumpairs,&riftspairs,numrifts,riftsnumsegments,riftssegments,x,y);
+
+	/*Order rifts so that they start from one tip, go to the other tip, and back: */
+	OrderRifts(&riftstips,riftssegments,riftspairs,numrifts,riftsnumsegments,x,y,nods,nel);
+
+	/*Create penalty pairs, used by Imp: */
+	PenaltyPairs(&riftspenaltypairs,&riftsnumpenaltypairs,numrifts,riftssegments,riftsnumsegments,riftspairs,riftstips,x,y);
+
+	/*Assign output pointers for mesh*/
+	*pnel            = nel;
+	*pindex          = index;
+	*px              = x;
+	*py              = y;
+	*pnods           = nods;
+	*psegments       = segments;
+	*psegmentmarkers = segmentmarkers;
+	*pnum_seg        = num_seg;
+
+	/*Assign output pointers for rifts*/
+	*pnumrifts             = numrifts;
+	*priftsnumsegments     = riftsnumsegments;
+	*priftssegments        = riftssegments;
+	*priftsnumpairs        = riftsnumpairs;
+	*priftspairs           = riftspairs;
+	*priftstips            = riftstips;
+	*priftspenaltypairs    = riftspenaltypairs;
+	*priftsnumpenaltypairs = riftsnumpenaltypairs;
+}