Index: /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp	(revision 7144)
+++ /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp	(revision 7145)
@@ -30,4 +30,12 @@
 	Vec    vec_incontour=NULL;
 	double*    incontour=NULL;
+
+	/*threading: */
+	InterpFromMesh2dxThreadStruct gate;
+	int num=1;
+
+	#ifdef _MULTITHREADING_
+	num=_NUMTHREADS_;
+	#endif
 
 	/*some checks*/
@@ -80,61 +88,26 @@
 	}
 
-	/*Loop over the elements*/
-	if (debug) printf("      interpolation progress:   %5.2lf %%",0.0);
-	for (i=0;i<nels_data;i++){
+	/*initialize thread parameters: */
+	gate.interpolation_type=interpolation_type;
+	gate.debug=debug;
+	gate.nels_data=nels_data;
+	gate.index_data=index_data;
+	gate.x_data=x_data;
+	gate.y_data=y_data;
+	gate.data=data;
+	gate.xmin=xmin;
+	gate.xmax=xmax;
+	gate.ymin=ymin;
+	gate.ymax=ymax;
+	gate.nods_prime=nods_prime;
+	gate.data_prime=data_prime;
+	gate.x_prime=x_prime;
+	gate.y_prime=y_prime;
+	gate.default_values=default_values;
+	gate.num_default_values=num_default_values;
+	gate.incontour=incontour;
 
-		/*display current iteration*/
-		if (debug && fmod((double)i,(double)100)==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)i/nels_data*100);
-
-		/*if there is no point inside the domain, go to next iteration*/
-		if ( (x_data[(int)index_data[3*i+0]-1]<xmin) && (x_data[(int)index_data[3*i+1]-1]<xmin) && (x_data[(int)index_data[3*i+2]-1]<xmin)) continue;
-		if ( (x_data[(int)index_data[3*i+0]-1]>xmax) && (x_data[(int)index_data[3*i+1]-1]>xmax) && (x_data[(int)index_data[3*i+2]-1]>xmax)) continue;
-		if ( (y_data[(int)index_data[3*i+0]-1]<ymin) && (y_data[(int)index_data[3*i+1]-1]<ymin) && (y_data[(int)index_data[3*i+2]-1]<ymin)) continue;
-		if ( (y_data[(int)index_data[3*i+0]-1]>ymax) && (y_data[(int)index_data[3*i+1]-1]>ymax) && (y_data[(int)index_data[3*i+2]-1]>ymax)) continue;
-
-		/*get area of the current element (Jacobian = 2 * area)*/
-		//area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
-		area=x_data[(int)index_data[3*i+1]-1]*y_data[(int)index_data[3*i+2]-1]-y_data[(int)index_data[3*i+1]-1]*x_data[(int)index_data[3*i+2]-1]
-		  +  x_data[(int)index_data[3*i+0]-1]*y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+0]-1]*x_data[(int)index_data[3*i+1]-1]
-		  +  x_data[(int)index_data[3*i+2]-1]*y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1]*x_data[(int)index_data[3*i+0]-1];
-
-		/*loop over the prime nodes*/
-		for (j=0;j<nods_prime;j++){
-
-			if(incontour[j]){
-
-				/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
-				area_1=((x_prime[j]-x_data[(int)index_data[3*i+2]-1])*(y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+2]-1]) 
-						-  (y_prime[j]-y_data[(int)index_data[3*i+2]-1])*(x_data[(int)index_data[3*i+1]-1]-x_data[(int)index_data[3*i+2]-1]))/area;
-				/*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
-				area_2=((x_data[(int)index_data[3*i+0]-1]-x_data[(int)index_data[3*i+2]-1])*(y_prime[j]-y_data[(int)index_data[3*i+2]-1]) 
-						- (y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1])*(x_prime[j]-x_data[(int)index_data[3*i+2]-1]))/area;
-				/*Get third area coordinate = 1-area1-area2*/
-				area_3=1-area_1-area_2;
-
-				/*is the current point in the current element?*/
-				if (area_1>=0 && area_2>=0 && area_3>=0){
-
-					/*Yes ! compute the value on the point*/
-					if (interpolation_type==1){
-						/*nodal interpolation*/
-						data_value=area_1*data[(int)index_data[3*i+0]-1]+area_2*data[(int)index_data[3*i+1]-1]+area_3*data[(int)index_data[3*i+2]-1];
-					}
-					else{
-						/*element interpolation*/
-						data_value=data[i];
-					}
-					if (isnan(data_value)){
-						if(num_default_values==1) data_value=default_values[0];
-						else data_value=default_values[j];
-					}
-
-					/*insert value and go to the next point*/
-					VecSetValue(data_prime,j,data_value,INSERT_VALUES);
-				}
-			}
-		}
-	}
-	 if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
+	/*launch the thread manager with InterpFromGridToMeshxt as a core: */
+	LaunchThread(InterpFromMesh2dxt,(void*)&gate,num);
 
 	/*Assign output pointers:*/
Index: /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h
===================================================================
--- /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h	(revision 7144)
+++ /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h	(revision 7145)
@@ -9,7 +9,33 @@
 #include "../../toolkits/toolkits.h"
 
+
+/*threading: */
+typedef struct{
+
+	int interpolation_type;
+	bool debug;
+	int  nels_data;
+	double* index_data;
+	double* x_data;
+	double* y_data;
+	double* data;
+	double xmin,xmax;
+	double ymin,ymax;
+	int    nods_prime;
+	Vec    data_prime;
+	double* x_prime;
+	double* y_prime;
+	double* default_values;
+	int     num_default_values;
+	double*    incontour;
+
+
+} InterpFromMesh2dxThreadStruct;
+
 int InterpFromMesh2dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime,
 		double* default_values,int num_default_values,Contour** contours,int numcontours);
 
+void* InterpFromMesh2dxt(void* vInterpFromMesh2dxThreadStruct);
+
 #endif /* _INTERPFROMMESH2DX_H */
 
Index: /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp
===================================================================
--- /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp	(revision 7145)
+++ /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp	(revision 7145)
@@ -0,0 +1,129 @@
+/*!\file:  InterpFromMesh2dxt.cpp
+ * \brief  thread core for InterpFromMesh2dxt code
+ */ 
+
+#include "./InterpFromMesh2dx.h"
+#include "../../shared/shared.h"
+
+void* InterpFromMesh2dxt(void* vpthread_handle){
+
+	/*gate variables :*/
+	InterpFromMesh2dxThreadStruct* gate=NULL;
+	pthread_handle* handle=NULL;
+	int     my_thread;
+	int     num_threads;
+	
+	int interpolation_type;
+	bool debug;
+	int  nels_data;
+	double* index_data=NULL;
+	double* x_data=NULL;
+	double* y_data=NULL;
+	double* data=NULL;
+	double xmin,xmax;
+	double ymin,ymax;
+	int    nods_prime;
+	Vec    data_prime=NULL;
+	double* x_prime=NULL;
+	double* y_prime=NULL;
+	double* default_values=NULL;
+	int     num_default_values;
+	double*    incontour=NULL;
+
+	/*intermediary: */
+	int     i0;
+	int     i1;
+	int     i,j;
+	double  area;
+	double  area_1,area_2,area_3;
+	double  data_value;
+
+	/*recover handle and gate: */
+	handle=(pthread_handle*)vpthread_handle;
+	gate=(InterpFromMesh2dxThreadStruct*)handle->gate;
+	my_thread=handle->id;
+	num_threads=handle->num;
+	
+	/*recover parameters :*/
+	interpolation_type=gate->interpolation_type;
+	debug=gate->debug;
+	nels_data=gate->nels_data;
+	index_data=gate->index_data;
+	x_data=gate->x_data;
+	y_data=gate->y_data;
+	data=gate->data;
+	xmin=gate->xmin;
+	xmax=gate->xmax;
+	ymin=gate->ymin;
+	ymax=gate->ymax;
+	nods_prime=gate->nods_prime;
+	data_prime=gate->data_prime;
+	x_prime=gate->x_prime;
+	y_prime=gate->y_prime;
+	default_values=gate->default_values;
+	num_default_values=gate->num_default_values;
+	incontour=gate->incontour;
+
+	/*partition loop across threads: */
+	PartitionRange(&i0,&i1,nels_data,num_threads,my_thread);
+
+	/*Loop over the elements*/
+	if (debug && my_thread==0) printf("      interpolation progress:   %5.2lf %%",0.0);
+
+	for (i=i0;i<i1;i++){
+
+		/*display current iteration*/
+		if (debug && my_thread==0 && fmod((double)i,(double)100)==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)i/nels_data*100*num_threads);
+
+		/*if there is no point inside the domain, go to next iteration*/
+		if ( (x_data[(int)index_data[3*i+0]-1]<xmin) && (x_data[(int)index_data[3*i+1]-1]<xmin) && (x_data[(int)index_data[3*i+2]-1]<xmin)) continue;
+		if ( (x_data[(int)index_data[3*i+0]-1]>xmax) && (x_data[(int)index_data[3*i+1]-1]>xmax) && (x_data[(int)index_data[3*i+2]-1]>xmax)) continue;
+		if ( (y_data[(int)index_data[3*i+0]-1]<ymin) && (y_data[(int)index_data[3*i+1]-1]<ymin) && (y_data[(int)index_data[3*i+2]-1]<ymin)) continue;
+		if ( (y_data[(int)index_data[3*i+0]-1]>ymax) && (y_data[(int)index_data[3*i+1]-1]>ymax) && (y_data[(int)index_data[3*i+2]-1]>ymax)) continue;
+
+		/*get area of the current element (Jacobian = 2 * area)*/
+		//area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
+		area=x_data[(int)index_data[3*i+1]-1]*y_data[(int)index_data[3*i+2]-1]-y_data[(int)index_data[3*i+1]-1]*x_data[(int)index_data[3*i+2]-1]
+		  +  x_data[(int)index_data[3*i+0]-1]*y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+0]-1]*x_data[(int)index_data[3*i+1]-1]
+		  +  x_data[(int)index_data[3*i+2]-1]*y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1]*x_data[(int)index_data[3*i+0]-1];
+
+		/*loop over the prime nodes*/
+		for (j=0;j<nods_prime;j++){
+
+			if(incontour[j]){
+
+				/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
+				area_1=((x_prime[j]-x_data[(int)index_data[3*i+2]-1])*(y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+2]-1]) 
+						-  (y_prime[j]-y_data[(int)index_data[3*i+2]-1])*(x_data[(int)index_data[3*i+1]-1]-x_data[(int)index_data[3*i+2]-1]))/area;
+				/*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
+				area_2=((x_data[(int)index_data[3*i+0]-1]-x_data[(int)index_data[3*i+2]-1])*(y_prime[j]-y_data[(int)index_data[3*i+2]-1]) 
+						- (y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1])*(x_prime[j]-x_data[(int)index_data[3*i+2]-1]))/area;
+				/*Get third area coordinate = 1-area1-area2*/
+				area_3=1-area_1-area_2;
+
+				/*is the current point in the current element?*/
+				if (area_1>=0 && area_2>=0 && area_3>=0){
+
+					/*Yes ! compute the value on the point*/
+					if (interpolation_type==1){
+						/*nodal interpolation*/
+						data_value=area_1*data[(int)index_data[3*i+0]-1]+area_2*data[(int)index_data[3*i+1]-1]+area_3*data[(int)index_data[3*i+2]-1];
+					}
+					else{
+						/*element interpolation*/
+						data_value=data[i];
+					}
+					if (isnan(data_value)){
+						if(num_default_values==1) data_value=default_values[0];
+						else data_value=default_values[j];
+					}
+
+					/*insert value and go to the next point*/
+					VecSetValue(data_prime,j,data_value,INSERT_VALUES);
+				}
+			}
+		}
+	}
+	if (debug && my_thread==0) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
+
+}