Index: /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp
===================================================================
--- /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp	(revision 5353)
+++ /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp	(revision 5354)
@@ -1,6 +1,6 @@
-/*!\file InterpFromMeshToMesh2dx
+/*!\file TriaSearchx
  */
 
-#include "./InterpFromMeshToMesh2dx.h"
+#include "./TriaSearchx.h"
 
 #include "../../shared/shared.h"
@@ -13,9 +13,8 @@
 using namespace std;
 
-int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
-			double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values, Contour** contours, int numcontours){
-	
+void TriaSearchx(double* ptria,double* index,int nel, double* x, double* y, int nods,double x0, double y0){
+
 	/*Output*/
-	double* data_interp=NULL;
+	double tria;
 
 	/*Intermediary*/
@@ -23,5 +22,4 @@
 	I2     I;
 	int    i,j,k;
-	int    it;
 	int    i0,i1,i2;
 	double areacoord[3];
@@ -31,116 +29,33 @@
 	int verbose=0;
 
-	/*default values: */
-	Vec    vec_incontour=NULL;
-	double*    incontour=NULL;
-	bool   skip_bamg=false;
-
-
-	/*Checks*/
-	if (data_cols<=0){
-		ISSMERROR("data provided has a negative number of columns");
-	}
-	if (data_rows!=nods_data && data_rows!=nels_data){
-		ISSMERROR("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows);
-	}
-	if((num_default_values) && (data_cols>1)){
-		ISSMERROR("data provided can only have 1 column if a default value is provided");
-	}
-	
-	/*If default values supplied, figure out which nodes are inside the contour, including the border of the contour: */
-	if(num_default_values){
-		ContourToNodesx( &vec_incontour,x_interp,y_interp,nods_interp,contours,numcontours,1);
-		VecToMPISerial(&incontour,vec_incontour);
-	}
-
-	/*Initialize output*/
-	if (verbose) printf("Initializing output vector\n");
-	data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double));
-
 	// read background mesh 
 	if (verbose) printf("Reading mesh\n");
-	Mesh Th(index_data,x_data,y_data,nods_data,nels_data); 
+	Mesh Th(index,x,y,nods,nel); 
 	Th.CreateSingleVertexToTriangleConnectivity();
 
-	//Loop over output nodes
-	if (verbose) printf("Loop over the nodes\n");
-	for(i=0;i<nods_interp;i++){
-		
-		/*reset skip_bamg: */
-		skip_bamg=false;
+	//Get current point coordinates
+	r.x=x0; r.y=y0;
+	I=Th.R2ToI2(r);
 
-		/*figure out if we should skip bamg logic: */
-		if(num_default_values){
-			if(!incontour[i]){
-				/*This node is not inside the contour. Skip Bamg logic and apply default value.: */
-				skip_bamg=true;
-			}
-		}
+	//Find triangle holding r/I
+	Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
 
-		if(skip_bamg==false){
-
-			//Get current point coordinates
-			r.x=x_interp[i]; r.y=y_interp[i];
-			I2 I=Th.R2ToI2(r);
-
-			//Find triangle holding r/I
-			Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
-
-			// internal point 
-			if (tb.det>0){ 
-				//Area coordinate
-				areacoord[0]= (double) dete[0]/ tb.det;
-				areacoord[1]= (double) dete[1] / tb.det;
-				areacoord[2]= (double) dete[2] / tb.det;
-				//3 vertices of the triangle
-				i0=Th.GetId(tb[0]);
-				i1=Th.GetId(tb[1]);
-				i2=Th.GetId(tb[2]);
-				//triangle number
-				it=Th.GetId(tb);
-			}
-			//external point
-			else {
-				//Get closest adjacent triangle (inside the mesh)
-				AdjacentTriangle ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
-				int k=ta;
-				Triangle &tc=*(Triangle*)ta;
-				//Area coordinate
-				areacoord[VerticesOfTriangularEdge[k][1]] = aa;
-				areacoord[VerticesOfTriangularEdge[k][0]] = bb;
-				areacoord[OppositeVertex[k]] = 1 - aa -bb;
-				//3 vertices of the triangle
-				i0=Th.GetId(tc[0]);
-				i1=Th.GetId(tc[1]);
-				i2=Th.GetId(tc[2]);
-				//triangle number
-				it=Th.GetId(tc);
-			}
-			
-			if (data_rows==nods_data){
-				for (j=0;j<data_cols;j++){
-					data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
-				}
-			}
-			else{
-				for (j=0;j<data_cols;j++){
-					if (it<0 || it>=nels_data){
-						ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data);
-					}
-					data_interp[i*data_cols+j]=data[data_cols*it+j];
-				}
-			}
-		}
-		else{
-			if(num_default_values==1) data_interp[i]=default_values[0];
-			else data_interp[i]=default_values[i];
-		}
+	// internal point 
+	if (tb.det>0){ 
+		//Area coordinate
+		areacoord[0]= (double) dete[0]/ tb.det;
+		areacoord[1]= (double) dete[1] / tb.det;
+		areacoord[2]= (double) dete[2] / tb.det;
+		//3 vertices of the triangle
+		i0=Th.GetId(tb[0]);
+		i1=Th.GetId(tb[1]);
+		i2=Th.GetId(tb[2]);
+		//triangle number
+		tria=(double)Th.GetId(tb);
 	}
-
+	//external point
+	else tria=NAN;
+	
 	/*Assign output pointers:*/
-	if (verbose) printf("Assigning output\n");
-	*pdata_interp=data_interp;
-
-	/*No error return*/
-	return 1;
+	*ptria=tria;
 }
Index: /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h
===================================================================
--- /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h	(revision 5353)
+++ /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h	(revision 5354)
@@ -1,14 +1,13 @@
-/*!\file:  InterpFromMeshToMesh2dx.h
+/*!\file:  TriaSearchx.h
  * \brief header file for Bamg module
  */ 
 
-#ifndef _INTERPFROMMESHTOMESH2DX_H
-#define _INTERPFROMMESHTOMESH2DX_H
+#ifndef _TRIASEARCHX_H
+#define _TRIASEARCHX_H
 
 #include "../../objects/objects.h"
 
 /* local prototypes: */
-int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
-			double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values,Contour** contours, int numcontours );
+void TriaSearchx(double* ptria,double* index,int nel, double* x, double* y, int nods,double x0, double y0);
 
 #endif
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 5353)
+++ /issm/trunk/src/mex/Makefile.am	(revision 5354)
@@ -60,4 +60,5 @@
 				SystemMatrices\
 				Test\
+				TriaSearch\
 				TriMesh\
 				TriMeshNoDensity\
@@ -260,4 +261,7 @@
 			  SystemMatrices/SystemMatrices.h
 
+TriaSearch_SOURCES = TriaSearch/TriaSearch.cpp\
+			  TriaSearch/TriaSearch.h
+
 TriMesh_SOURCES = TriMesh/TriMesh.cpp\
 			  TriMesh/TriMesh.h
Index: /issm/trunk/src/mex/TriaSearch/TriaSearch.cpp
===================================================================
--- /issm/trunk/src/mex/TriaSearch/TriaSearch.cpp	(revision 5353)
+++ /issm/trunk/src/mex/TriaSearch/TriaSearch.cpp	(revision 5354)
@@ -1,6 +1,6 @@
-/*\file InterpFromMeshToMesh2d.c
- *\brief: bamg module.
+/*\file TriaSearch.c
+ *\brief: TriaSearch module. See TriaSearchx for more details.
  */
-#include "./InterpFromMeshToMesh2d.h"
+#include "./TriaSearch.h"
 
 void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
@@ -10,39 +10,16 @@
 	/*input: */
 	double* index=NULL;
-	int     index_cols;
+	int     nel;
+	int     dummy;
 
-	double* x_data=NULL;
-	int     x_data_rows;
+	double* x=NULL;
+	int     nods;
 
-	double* y_data=NULL;
-	int     y_data_rows;
+	double* y=NULL;
 
-	double* data=NULL; 
-	int     data_rows;
-	int     data_cols;
-
-	double* x_interp=NULL;
-	double* y_interp=NULL;
-
-	int     x_interp_rows;
-	int     y_interp_rows;
-
-	double* default_values=NULL;
-	int     num_default_values=0;
-
-	//contours
-	mxArray*  matlabstructure=NULL;
-	Contour** contours=NULL;
-	int       numcontours;
-	Contour*  contouri=NULL;
-
-	/*Intermediary*/
-	int nods_data;
-	int nels_data;
-	int nods_interp;
-	int verbose=0;
-
+	double  x0,y0;
+	
 	/* output: */
-	double* data_interp=NULL;
+	double  tria;
 
 	/*Boot module: */
@@ -50,84 +27,25 @@
 
 	/*checks on arguments on the matlab side: */
-	if(nlhs!=NLHS){
-		InterpFromMeshToMesh2dUsage();
-		ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
-	}
-	if((nrhs!=6) & (nrhs!=8)){
-		InterpFromMeshToMesh2dUsage();
-		ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
-	}
+	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&TriaSearchUsage);
 
 	/*Input datasets: */
-	if (verbose) printf("Fetching inputs\n");
-	FetchData(&index,&nels_data,&index_cols,INDEXHANDLE);
-	FetchData(&x_data,&x_data_rows,NULL,XHANDLE);
-	FetchData(&y_data,&y_data_rows,NULL,YHANDLE);
-	FetchData(&data,&data_rows,&data_cols,DATAHANDLE);
-	FetchData(&x_interp,&x_interp_rows,NULL,XINTERPHANDLE);
-	FetchData(&y_interp,&y_interp_rows,NULL,YINTERPHANDLE);
+	FetchData(&index,&nel,&dummy,INDEXHANDLE);
+	FetchData(&x,&nods,XHANDLE);
+	FetchData(&y,&nods,YHANDLE);
+	FetchData(&x0,X0HANDLE);
+	FetchData(&y0,Y0HANDLE);
 
-	if(nrhs==8){
-		
-		/*Call expread on filename to build a contour array in the matlab workspace: */
-		mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread");
-
-		/*default values: */
-		FetchData(&default_values,&num_default_values,DEFAULTHANDLE);
-
-		/*contours: */
-		numcontours=mxGetNumberOfElements(matlabstructure);
-		contours=(Contour**)xmalloc(numcontours*sizeof(Contour*));
-		for(i=0;i<numcontours;i++){
-			//allocate
-			contouri=(Contour*)xmalloc(sizeof(Contour));
-			//retrieve dimension of this contour.
-			contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
-			//set pointers.
-			contouri->x=mxGetPr(mxGetField(matlabstructure,i,"x"));
-			contouri->y=mxGetPr(mxGetField(matlabstructure,i,"y"));
-			*(contours+i)=contouri;
-		}
-
-		/* Debugging of contours :{{{1*/
-		/*for(i=0;i<numcontours;i++){
-		  printf("\nContour echo: contour number  %i / %i\n",i+1,numcontours);
-		  contouri=*(contours+i);
-		  printf("   Number of grids %i\n",contouri->nods);
-		  for (j=0;j<contouri->nods;j++){
-		  printf("   %lf %lf\n",*(contouri->x+j),*(contouri->y+j));
-		  }
-		  }*/
-		/*}}}*/
-	}
-	else{
-		num_default_values=0;
-		default_values=NULL;
-		numcontours=0;
-		contours=NULL;
-	}
-
-	/*some checks*/
-	if (verbose) printf("Checking inputs\n");
-	if (x_data_rows!=y_data_rows){
-		ISSMERROR("vectors x and y should have the same length!");
-	}
-	if (x_interp_rows!=y_interp_rows){
-		ISSMERROR("vectors x_interp and y_interp should have the same length!");
-	}
-	if (index_cols!=3){
-		ISSMERROR("index should have 3 columns (input provided has %i columns)",index_cols);
-	}
-
-	/*get number of elements and number of nodes in the data*/
-	nods_data=x_data_rows;
-	nods_interp=x_interp_rows;
+	/* Echo: {{{1*/
+	//printf("(x0,y0)=(%g,%g)\n",x0,y0);
+	/*}}}*/
 
 	/* Run core computations: */
-	if (verbose) printf("Call core\n");
-	InterpFromMeshToMesh2dx(&data_interp,index,x_data,y_data,nods_data,nels_data,data,data_rows,data_cols,x_interp,y_interp,nods_interp,default_values,num_default_values,contours,numcontours);
+	TriaSearchx(&tria,index,nel,x,y,nods,x0,y0);
+
+	/* c to matlab: */
+	tria++;
 
 	/*Write data: */
-	WriteData(DATAINTERP,data_interp,nods_interp,data_cols);
+	WriteData(TRIA,tria);
 
 	/*end module: */
@@ -135,28 +53,12 @@
 }
 
-void InterpFromMeshToMesh2dUsage(void)
+void TriaSearchUsage(void)
 {
-	_printf_("INTERFROMMESHTOMESH2D - interpolation from a 2d triangular mesh onto a list of point\n");
-	_printf_("\n");
-	_printf_("   This function is a multi-threaded mex file that interpolates a field\n");
-	_printf_("   defined on a triangular mesh onto a list of point\n");
+	_printf_("TriaSearch- find triangle holding a point (x0,y0) in a mesh\n");
 	_printf_("\n");
 	_printf_("   Usage:\n");
-	_printf_("         data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp);\n");
-	_printf_("      or data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value,contourname);\n");
-	_printf_("\n");
-	_printf_("      index: index of the mesh where data is defined\n");
-	_printf_("      x,y: coordinates of the nodes where data is defined\n");
-	_printf_("      data: matrix holding the data to be interpolated onto the mesh. (one column per field)\n");
-	_printf_("      x_interp,y_interp: coordinates of the points onto which we interpolate.\n");
-	_printf_("      if default_value and contourname not specified: linear interpolation will happen on all x_interp,y_interp.\n");
-	_printf_("      if (default_value,contourname) specified: linear interpolation will happen on all x_interp,y_interp inside the contour, default value will be adopted on the rest of the mesh.\n");
-	_printf_("      note that default_value is either a scalar, or a vector of size  length(x_interp)\n");
-	_printf_("      data_interp: vector of mesh interpolated data.\n");
-	_printf_("\n");
-	_printf_("   Example:\n");
-	_printf_("      load('temperature.mat');\n");
-	_printf_("      md.temperature=InterpFromMeshToMesh2d(index,x,y,temperature,md.x,md.y);\n");
-	_printf_("      md.temperature=InterpFromMeshToMesh2d(index,x,y,temperature,md.x,md.y,253,'Contour.exp');\n");
+	_printf_("         tria=TriaSearch(index,x,y,x0,y0);\n");
+	_printf_("      index,x,y: mesh triangulatrion\n");
+	_printf_("      x0,y0: coordinates of the point for which we are trying to find a triangle\n");
 	_printf_("\n");
 }
Index: /issm/trunk/src/mex/TriaSearch/TriaSearch.h
===================================================================
--- /issm/trunk/src/mex/TriaSearch/TriaSearch.h	(revision 5353)
+++ /issm/trunk/src/mex/TriaSearch/TriaSearch.h	(revision 5354)
@@ -1,11 +1,10 @@
-/*!\file InterpFromMeshToMesh2d.h
- * \brief: prototype for Data Interpolation mex module.
+/*!\file TriaSearch.h
  */
 
-#ifndef _INTERPFROMMESHTOMESH2d_H
-#define _INTERPFROMMESHTOMESH2d_H
+#ifndef _TRIASEARCH_H
+#define _TRIASEARCH_H
 
 /* local prototypes: */
-void InterpFromMeshToMesh2dUsage(void);
+void TriaSearchUsage(void);
 
 #include "../../c/modules/modules.h"
@@ -14,5 +13,5 @@
 
 #undef __FUNCT__ 
-#define __FUNCT__  "InterpFromMeshToMesh2d"
+#define __FUNCT__  "TriaSearch"
 
 
@@ -21,12 +20,9 @@
 #define XHANDLE prhs[1]
 #define YHANDLE prhs[2]
-#define DATAHANDLE prhs[3]
-#define XINTERPHANDLE prhs[4]
-#define YINTERPHANDLE prhs[5]
-#define DEFAULTHANDLE prhs[6]
-#define FILENAME prhs[7]
+#define X0HANDLE prhs[3]
+#define Y0HANDLE prhs[4]
 
 /* serial output macros: */
-#define DATAINTERP (mxArray**)&plhs[0]
+#define TRIA (mxArray**)&plhs[0]
 
 /* serial arg counts: */
@@ -34,3 +30,6 @@
 #define NLHS  1
 
+#undef NRHS
+#define NRHS  5
+
 #endif