Index: /issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridx.cpp
===================================================================
--- /issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridx.cpp	(revision 2285)
+++ /issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridx.cpp	(revision 2285)
@@ -0,0 +1,167 @@
+/*!\file:  InterpFromGridx.cpp
+ * \brief  "c" core code for interpolating values from a structured grid.
+ */ 
+
+#include "./InterpFromGridx.h"
+#include "../shared/shared.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "InterpFromGridx"
+
+int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
+
+int InterpFromGridx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value) {
+
+
+	/*output: */
+	Vec data_mesh=NULL;
+	
+	/*Intermediary*/
+	double* x=NULL;
+	double* y=NULL;
+	int i,m,n;
+	double x_grid,y_grid;
+	double xi,eta;
+	double G1,G2,G3,G4,data_value;
+	double area_1,area_2,area_3;
+	double area;
+
+	/*Some checks on arguments: */
+	if ((M<=2) || (N<=2) || (nods<=0)){
+		throw ErrorException(__FUNCT__,"nothing to be done according to the dimensions of input matrices and vectors.");
+	}
+
+	/*Allocate output vector: */
+	data_mesh=NewVec(nods);
+
+	/*Find out what kind of coordinates (x_in,y_in) have been given is input*/
+	if(N==(x_rows-1) && M==(y_rows-1)){
+
+		/*The coordinates given in input describe the contour of each pixel. Take the center of each pixel*/
+		x=(double*)xmalloc(N*sizeof(double));
+		y=(double*)xmalloc(M*sizeof(double));
+		for (i=0;i<N;i++) x[i]=(x_in[i]+x_in[i+1])/2;
+		for (i=0;i<M;i++) y[i]=(y_in[i]+y_in[i+1])/2;
+		x_rows=x_rows-1;
+		y_rows=y_rows-1;
+	}
+	else if (M==x_rows && N==x_rows){
+
+		/*The coordinates given in input describe the center each pixel. Keep them*/
+		x=(double*)xmalloc(N*sizeof(double));
+		y=(double*)xmalloc(M*sizeof(double));
+		for (i=0;i<N;i++) x[i]=x_in[i];
+		for (i=0;i<M;i++) y[i]=y_in[i];
+	}
+	else{
+		throw ErrorException(__FUNCT__,"x and y vectors length should be 1 or 0 more than data number of rows.");
+	}
+
+	/*Linear (triangle) interpolation: */
+	for ( i=MPI_Lowerrow(nods); i<MPI_Upperrow(nods); i++) {
+
+		x_grid=*(x_mesh+i);
+		y_grid=*(y_mesh+i);
+
+		/*Find indices m and n into y and x, for which  y(m)<=y_grids<=y(m+1) and x(n)<=x_grid<=x(n+1)*/
+		if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
+			
+			/*Get area*/
+			area=(x[n+1]-x[n])*(y[m+1]-y[m]);
+
+			/*is it the upper right triangle?*/
+			/*2'     3'
+			 *+-----+
+			 *1\    |
+			 *| \   |
+			 *|  \  |
+			 *|   \ |
+			 *|    \|
+			 *2----3+1' */
+			if ((x_grid-x[n])/(x[n+1]-x[n])<(y_grid-y[m])/(y[m+1]-y[m])){
+
+				/*Then find values of data at each summit*/
+				G1=*(data+m*N+n);
+				G2=*(data+(m+1)*N+n+1);
+				G3=*(data+(m+1)*N+n);
+
+				/*Get first area coordinate*/
+				area_1=((y[m+1]-y_grid)*(x[n+1]-x[n]))/area;
+				/*Get second area coordinate*/
+				area_2=((x_grid-x[n])*(y[m+1]-y[m]))/area;
+				/*Get third area coordinate = 1-area1-area2*/
+				area_3=1-area_1-area_2;
+
+				/*interpolate*/
+				data_value=area_1*G1+area_2*G2+area_3*G3;
+			}
+			else {
+
+				/*Then find values of data at each summit*/
+				G1=*(data+(m+1)*N+n+1);
+				G2=*(data+m*N+n);
+				G3=*(data+m*N+n+1);
+
+				/*Get first area coordinate*/
+				area_1=((y_grid-y[m])*(x[n+1]-x[n]))/area;
+				/*Get second area coordinate*/
+				area_2=((x[n+1]-x_grid)*(y[m+1]-y[m]))/area;
+				/*Get third area coordinate = 1-area1-area2*/
+				area_3=1-area_1-area_2;
+
+				/*interpolate*/
+				data_value=area_1*G1+area_2*G2+area_3*G3;
+			}
+
+			/*Treat NANs*/
+			if (isnan(data_value)){
+				data_value=default_value;
+			}
+		}
+		else{
+			data_value=default_value;
+		}
+		VecSetValue(data_mesh,i,data_value,INSERT_VALUES);
+	}
+
+	/*Assign output pointers:*/
+	*pdata_mesh=data_mesh;
+}
+
+int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid){
+
+	int foundx=0;
+	int foundy=0;
+	int i;
+	int m=-1;
+	int n=-1;
+
+	for (i=0;i<x_rows-1;i++){
+		if ( (*(x+i)<=xgrid) && (xgrid<*(x+i+1)) ){
+			m=i;
+			foundx= 1;
+			break;
+		}
+	}
+	if(*(x+x_rows-1)==xgrid){
+		m=x_rows-2;
+		foundx=1;
+	}
+	
+	for (i=0;i<y_rows-1;i++){
+		if ( (*(y+i)<=ygrid) && (ygrid<*(y+i+1)) ){
+			n=i;
+			foundy= 1;
+			break;
+		}
+	}
+	if(*(y+y_rows-1)==ygrid){
+		m=y_rows-2;
+		foundy=1;
+	}
+
+	/*Assign output pointers:*/
+	*pm=m;
+	*pn=n;
+	return foundx*foundy;
+}
Index: /issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridx.h
===================================================================
--- /issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridx.h	(revision 2285)
+++ /issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridx.h	(revision 2285)
@@ -0,0 +1,13 @@
+/*!\file InterpFromGridx.h
+ * \brief: header file for Data interpolation routines.
+ */
+
+#ifndef _INTERPFROMGRIDX_H
+#define _INTERPFROMGRIDX_H
+
+#include "../toolkits/toolkits.h"
+
+int InterpFromGridx( Vec* pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value);
+
+#endif /* _INTERPFROMGRIDX_H */
+
Index: /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dx.cpp
===================================================================
--- /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dx.cpp	(revision 2285)
+++ /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dx.cpp	(revision 2285)
@@ -0,0 +1,112 @@
+/*!\file:  InterpFromMesh2dx.cpp
+ * \brief  "c" core code for interpolating values from a structured grid.
+ */ 
+
+#include "./InterpFromMesh2dx.h"
+#include "../shared/shared.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "InterpFromMesh2dx"
+
+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_value) {
+
+	/*Output*/
+	Vec data_prime=NULL;
+
+	/*Intermediary*/
+	int i,j;
+	int interpolation_type;
+	bool debug;
+	double area;
+	double area_1,area_2,area_3;
+	double data_value;
+	double xmin,xmax;
+	double ymin,ymax;
+
+	/*some checks*/
+	if (nels_data<1 || nods_data<3 || nods_prime==0){
+		throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
+	}
+
+	/*Set debug to 1 if there are lots of elements*/
+	debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9));
+
+	/*figure out what kind of interpolation is needed*/
+	if (data_length==nods_data){
+		interpolation_type=1;
+	}
+	else if (data_length==nels_data){
+		interpolation_type=2;
+	}
+	else{
+		throw ErrorException(__FUNCT__,"length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
+	}
+
+	/*Get prime mesh extrema coordinates*/
+	xmin=x_prime[0]; xmax=x_prime[0];ymin=y_prime[0]; ymax=y_prime[0];
+	for (i=1;i<nods_prime;i++){
+		if (x_prime[i]<xmin) xmin=x_prime[i];
+		if (x_prime[i]>xmax) xmax=x_prime[i];
+		if (y_prime[i]<ymin) ymin=y_prime[i];
+		if (y_prime[i]>ymax) ymax=y_prime[i];
+	}
+
+	/*Initialize output*/
+	data_prime=NewVec(nods_prime);
+	for (i=0;i<nods_prime;i++) VecSetValue(data_prime,i,default_value,INSERT_VALUES);
+
+	/*Loop over the elements*/
+	if (debug) printf("      interpolation progress:   %5.2lf %%",0.0);
+	for (i=0;i<nels_data;i++){
+
+		/*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++){
+
+			/*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)) data_value=default_value;
+
+				/*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);
+
+	/*Assign output pointers:*/
+	*pdata_prime=data_prime;
+}
Index: /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dx.h
===================================================================
--- /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dx.h	(revision 2285)
+++ /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dx.h	(revision 2285)
@@ -0,0 +1,13 @@
+/*!\file InterpFromMesh2dx.h
+ * \brief: header file for Data interpolation routines.
+ */
+
+#ifndef _INTERPFROMMESH2DX_H
+#define _INTERPFROMMESH2DX_H
+
+#include "../toolkits/toolkits.h"
+
+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_value);
+
+#endif /* _INTERPFROMMESH2DX_H */
+
Index: /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dx.cpp
===================================================================
--- /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dx.cpp	(revision 2285)
+++ /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dx.cpp	(revision 2285)
@@ -0,0 +1,122 @@
+/*!\file:  InterpFromMesh3dx.cpp
+ * \brief  "c" core code for interpolating values from a structured grid.
+ */ 
+
+#include "./InterpFromMesh3dx.h"
+#include "../shared/shared.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "InterpFromMesh3dx"
+
+int InterpFromMesh3dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) {
+
+	/*Output*/
+	Vec data_prime=NULL;
+
+	/*Intermediary*/
+	int i,j;
+	int interpolation_type;
+	bool debug;
+	double area;
+	double area_1,area_2,area_3;
+	double zeta,bed,surface;
+	double data_value;
+	double xmin,xmax;
+	double ymin,ymax;
+
+	/*some checks*/
+	if (nels_data<1 || nods_data<6 || nods_prime==0){
+		throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
+	}
+
+	/*Set debug to 1 if there are lots of elements*/
+	debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9));
+
+	/*figure out what kind of interpolation is needed*/
+	if (data_length==nods_data){
+		interpolation_type=1;
+	}
+	else if (data_length==nels_data){
+		interpolation_type=2;
+	}
+	else{
+		throw ErrorException(__FUNCT__,"length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
+	}
+
+	/*Get prime mesh extrema coordinates*/
+	xmin=x_prime[0]; xmax=x_prime[0];ymin=y_prime[0]; ymax=y_prime[0];
+	for (i=1;i<nods_prime;i++){
+		if (x_prime[i]<xmin) xmin=x_prime[i];
+		if (x_prime[i]>xmax) xmax=x_prime[i];
+		if (y_prime[i]<ymin) ymin=y_prime[i];
+		if (y_prime[i]>ymax) ymax=y_prime[i];
+	}
+
+	/*Initialize output*/
+	data_prime=NewVec(nods_prime);
+	for (i=0;i<nods_prime;i++) VecSetValue(data_prime,i,default_value,INSERT_VALUES);
+
+	/*Loop over the elements*/
+	if (debug) printf("      interpolation progress:   %5.2lf %%",0.0);
+	for (i=0;i<nels_data;i++){
+
+		/*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[6*i+0]-1]<xmin) && (x_data[(int)index_data[6*i+1]-1]<xmin) && (x_data[(int)index_data[6*i+2]-1]<xmin)) continue;
+		if ( (x_data[(int)index_data[6*i+0]-1]>xmax) && (x_data[(int)index_data[6*i+1]-1]>xmax) && (x_data[(int)index_data[6*i+2]-1]>xmax)) continue;
+		if ( (y_data[(int)index_data[6*i+0]-1]<ymin) && (y_data[(int)index_data[6*i+1]-1]<ymin) && (y_data[(int)index_data[6*i+2]-1]<ymin)) continue;
+		if ( (y_data[(int)index_data[6*i+0]-1]>ymax) && (y_data[(int)index_data[6*i+1]-1]>ymax) && (y_data[(int)index_data[6*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[6*i+1]-1]*y_data[(int)index_data[6*i+2]-1]-y_data[(int)index_data[6*i+1]-1]*x_data[(int)index_data[6*i+2]-1]
+		  +  x_data[(int)index_data[6*i+0]-1]*y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+0]-1]*x_data[(int)index_data[6*i+1]-1]
+		  +  x_data[(int)index_data[6*i+2]-1]*y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1]*x_data[(int)index_data[6*i+0]-1];
+
+		/*loop over the prime nodes*/
+		for (j=0;j<nods_prime;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[6*i+2]-1])*(y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+2]-1]) 
+						-  (y_prime[j]-y_data[(int)index_data[6*i+2]-1])*(x_data[(int)index_data[6*i+1]-1]-x_data[(int)index_data[6*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[6*i+0]-1]-x_data[(int)index_data[6*i+2]-1])*(y_prime[j]-y_data[(int)index_data[6*i+2]-1]) 
+						- (y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1])*(x_prime[j]-x_data[(int)index_data[6*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 2d element?*/
+			if (area_1>=0 && area_2>=0 && area_3>=0){
+
+				/*compute bottom and top height of the element at this 2d position*/
+				bed    =area_1*z_data[(int)index_data[6*i+0]-1]+area_2*z_data[(int)index_data[6*i+1]-1]+area_3*z_data[(int)index_data[6*i+2]-1];
+				surface=area_1*z_data[(int)index_data[6*i+3]-1]+area_2*z_data[(int)index_data[6*i+4]-1]+area_3*z_data[(int)index_data[6*i+5]-1];
+
+				/*Compute zeta*/
+				zeta=2*(z_prime[j]-bed)/(surface-bed)-1;
+
+				if (zeta >=-1 && zeta<=1){
+					if (interpolation_type==1){
+						/*nodal interpolation*/
+						data_value=(1-zeta)/2*(area_1*data[(int)index_data[6*i+0]-1]+area_2*data[(int)index_data[6*i+1]-1]+area_3*data[(int)index_data[6*i+2]-1]) 
+						  +        (1+zeta)/2*(area_1*data[(int)index_data[6*i+3]-1]+area_2*data[(int)index_data[6*i+4]-1]+area_3*data[(int)index_data[6*i+5]-1]);
+					}
+					else{
+						/*element interpolation*/
+						data_value=data[i];
+					}
+					if (isnan(data_value)) data_value=default_value;
+
+					/*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);
+
+	/*Assign output pointers:*/
+	*pdata_prime=data_prime;
+}
Index: /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dx.h
===================================================================
--- /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dx.h	(revision 2285)
+++ /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dx.h	(revision 2285)
@@ -0,0 +1,13 @@
+/*!\file InterpFromMesh3dx.h
+ * \brief: header file for Data interpolation routines.
+ */
+
+#ifndef _INTERPFROMMESH3DX_H
+#define _INTERPFROMMESH3DX_H
+
+#include "../toolkits/toolkits.h"
+
+int InterpFromMesh3dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value);
+
+#endif /* _INTERPFROMMESH3DX_H */
+