/*! \file  MeshProfileIntersection.cpp
    \brief: takes a  .exp file (made of several profiles), and figures out its intersection 
	with a mesh.

	usage:
	[segments]=MeshProfileIntersection(index,x,y,filename);

	where:
	input:
		index,x,y is a triangulation
		filename: name of Argus style .exp file containing the segments (can be groups of disconnected segments)
	output:
		segments: array made of x1,y1,x2,y2,element_id lines (x1,y1) and (x2,y2) are segment extremities for a segment 
		belonging to the elemnt_id element. there are as many lines in segments as there are segments intersecting the 
		mesh.
*/
	
#include "./MeshProfileIntersection.h"

void MeshProfileIntersectionUsage(void){/*{{{*/
	_printLine_("   usage:");
	_printLine_("   [segments]=MeshProfileIntersection(index,x,y,filename);");
	_printLine_("   where:");
	_printLine_("   input:");
	_printLine_("        index,x,y is a triangulation");
	_printLine_("        filename: name of Argus style .exp file containing the segments (can be groups of disconnected segments)");
	_printLine_("   output:");
	_printLine_("        segments: array made of x1,y1,x2,y2,element_id lines (x1,y1) and (x2,y2) are segment extremities for a segment ");
	_printLine_("        belonging to the elemnt_id element. there are as many lines in segments as there are segments intersecting the ");
	_printLine_("        mesh.");
}/*}}}*/
WRAPPER(MeshProfileIntersection){

	int i,j;

	/* required input: */
	//mesh
	double* double_index=NULL;
	int*    index=NULL;
	int     nel;
	double* x=NULL;
	double* y=NULL;
	int     nods;
	int     dummy;
	
	//contours
	DataSet *domain = NULL;
	Contour<double>** contours=NULL;
	int       numcontours;
	Contour<double>*  contouri=NULL;

	/* output: */
	double* segments=NULL;
	int     numsegs;

	/*Boot module: */
	MODULEBOOT();

	/*checks on arguments: */
	CHECKARGUMENTS(NLHS,NRHS,&MeshProfileIntersectionUsage);

	/*Fetch inputs: */
	//index
	FetchData(&double_index,&nel,&dummy,INDEX);
	if(dummy!=3 && dummy!=6)_error_("element triangulation should be of 3 or 6 column width!");
	index=xNew<int>(nel*3);
	for(i=0;i<nel;i++){
		for(j=0;j<3;j++){
			*(index+3*i+j)=(int)*(double_index+dummy*i+j)-1; //"C" style indexing
		}
	}
	//x and y
	FetchData(&x,&nods,X);
	FetchData(&y,&dummy,Y);

	//contours
	FetchData(&domain,FILENAME);
	// MeshProfileIntersectionx should be modified to take DataSet directly (and perhaps SeqMat and SeqVec).
	numcontours=domain->Size();
	contours=xNew<Contour<double>*>(numcontours);
	for(i=0;i<numcontours;i++)
		*(contours+i)=(Contour<double>*)domain->GetObjectByOffset(i);

	/* Debugging of contours :{{{1*/
	/*for(i=0;i<numcontours;i++){
		_printLine_("\nContour echo: contour number  " << i+1 << " / " << numcontours);
		contouri=*(contours+i);
		_printLine_("   Number of vertices " << contouri->nods);
		for (j=0;j<contouri->nods;j++){
			_printLine_("   " << *(contouri->x+j) << "f " << *(contouri->y+j) << "f");
		}
	}*/
	/*}}}*/

	/*Run interpolation routine: */
	MeshProfileIntersectionx(&segments,&numsegs,index,x,y,nel,nods,contours,numcontours);

	/* output: */
	WriteData(SEGMENTS,segments,numsegs,5);

	/*end module: */
	MODULEEND();

}
