Index: issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp
===================================================================
--- issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp	(revision 4981)
+++ issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp	(revision 5032)
@@ -5,4 +5,6 @@
 
 void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+
+	int i;
 
 	/*input: */
@@ -26,5 +28,12 @@
 	int     y_interp_rows;
 
-	double default_value;
+	double* default_values=NULL;
+	int     num_default_values=0;
+
+	//contours
+	mxArray*  matlabstructure=NULL;
+	Contour** contours=NULL;
+	int       numcontours;
+	Contour*  contouri=NULL;
 
 	/*Intermediary*/
@@ -41,5 +50,12 @@
 
 	/*checks on arguments on the matlab side: */
-	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InterpFromMeshToMesh2dUsage);
+	if(nlhs!=NLHS){
+		InterpFromMeshToMesh2dUsage();
+		ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
+	}
+	if((nrhs!=6) & (nrhs!=8)){
+		InterpFromMeshToMesh2dUsage();
+		ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
+	}
 
 	/*Input datasets: */
@@ -51,5 +67,44 @@
 	FetchData(&x_interp,&x_interp_rows,NULL,XINTERPHANDLE);
 	FetchData(&y_interp,&y_interp_rows,NULL,YINTERPHANDLE);
-	FetchData(&default_value,DEFAULTHANDLE);
+
+	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*/
@@ -71,5 +126,5 @@
 	/* 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_value);
+	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);
 
 	/*Write data: */
@@ -88,5 +143,6 @@
 	_printf_("\n");
 	_printf_("   Usage:\n");
-	_printf_("      data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value);\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");
@@ -94,5 +150,7 @@
 	_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_("      default_value: default value if no data is found (holes).\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");
