Index: /issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp
===================================================================
--- /issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp	(revision 2284)
+++ /issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp	(revision 2284)
@@ -0,0 +1,164 @@
+/*!\file:  InterpFromMeshToGridx.cpp
+ * \brief  "c" core code for interpolating values from a structured grid.
+ */ 
+
+#include "./InterpFromMeshToGridx.h"
+#include "../shared/shared.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "InterpFromMeshToGridx"
+
+int InterpFromMeshToGridx(double** px_m,double** py_m,double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh, int data_length, double cornereast,double cornernorth,double xposting,double yposting,int nlines,int ncols,double default_value) {
+
+	/*Output*/
+	double* griddata=NULL;
+	double* x_m=NULL;
+	double* y_m=NULL;
+
+	/*Intermediary*/
+	int    i,j,n;
+	int    interpolation_type;
+	bool   debug;
+	double area;
+	double area_1,area_2,area_3;
+	double x_tria_min,y_tria_min;
+	double x_tria_max,y_tria_max;
+	double x_grid_min,y_grid_min;
+	double x_grid_max,y_grid_max;
+	double data_value;
+	double cornersouth,cornerwest;
+	double* x_grid=NULL;
+	double* y_grid=NULL;
+
+	/*some checks*/
+	if (nels<1 || nods<3 || nlines<1 || ncols<1 || xposting==0 || yposting==0){
+		throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
+	}
+
+	/*figure out what kind of interpolation is needed*/
+	if (data_length==nods){
+		interpolation_type=1;
+	}
+	else if (data_length==nels){
+		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)!");
+	}
+
+	/*First, allocate pointers: */
+	griddata=(double*)xcalloc(nlines*ncols,sizeof(double));
+	x_m=(double*)xcalloc(ncols+1,sizeof(double));
+	y_m=(double*)xcalloc(nlines+1,sizeof(double));
+	x_grid=(double*)xcalloc(ncols,sizeof(double));
+	y_grid=(double*)xcalloc(nlines,sizeof(double));
+
+	/*Set debug to 1 if there are lots of elements*/
+	debug=(bool)((double)nels >= pow((double)10,(double)4));
+
+	/*Initialize coordintes and griddata*/
+	// EAST = ACTUAL WEST !!!!!!!!!!!!!!!
+	cornersouth=cornernorth-nlines*yposting;
+	cornerwest =cornereast +ncols *xposting;
+	for(i=0;i<nlines;i++) y_grid[i]= cornersouth + yposting/2 + yposting*i;
+	for(i=0;i<ncols; i++) x_grid[i]= cornereast  + xposting/2 + xposting*i;
+	for(i=0;i<=nlines;i++)   y_m[i]   = cornersouth + yposting*i;
+	for(i=0;i<=ncols; i++)   x_m[i]   = cornereast  + xposting*i;
+	for(i=0;i<nlines;i++){
+		for(j=0;j<ncols; j++){
+			griddata[i*ncols+j]=default_value;
+		}
+	}
+	
+	/*Get extreme coordinates of the grid*/
+	if (xposting>0){
+		x_grid_min=cornereast;
+		x_grid_max=cornerwest;
+	}
+	else{
+		x_grid_min=cornerwest;
+		x_grid_max=cornereast;
+	}
+	if (yposting>0){
+		y_grid_min=cornersouth;
+		y_grid_max=cornernorth;
+	}
+	else{
+		y_grid_min=cornernorth;
+		y_grid_max=cornersouth;
+	}
+
+	/*Loop over the elements*/
+	if (debug) printf("      interpolation progress:   %5.2lf %%",0.0);
+	for (n=0;n<nels;n++){
+
+		/*display current iteration*/
+		if (debug && fmod((double)n,(double)100)==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)i/nels*100);
+
+		/*Get extrema coordinates of current elements*/
+		x_tria_min=x_mesh[(int)index_mesh[3*n+0]-1]; x_tria_max=x_tria_min;
+		y_tria_min=y_mesh[(int)index_mesh[3*n+0]-1]; y_tria_max=y_tria_min;
+		for (i=1;i<3;i++){
+			if(x_mesh[(int)index_mesh[3*n+i]-1]<x_tria_min) x_tria_min=x_mesh[(int)index_mesh[3*n+i]-1];
+			if(x_mesh[(int)index_mesh[3*n+i]-1]>x_tria_max) x_tria_max=x_mesh[(int)index_mesh[3*n+i]-1];
+			if(y_mesh[(int)index_mesh[3*n+i]-1]<y_tria_min) y_tria_min=y_mesh[(int)index_mesh[3*n+i]-1];
+			if(y_mesh[(int)index_mesh[3*n+i]-1]>y_tria_max) y_tria_max=y_mesh[(int)index_mesh[3*n+i]-1];
+		}
+
+		/*if the current triangle is not in the grid, continue*/
+		if ( (x_tria_min>x_grid_max) || (x_tria_max<x_grid_min) || (y_tria_min>y_grid_max) || (y_tria_max<y_grid_min) ) 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_mesh[(int)index_mesh[3*n+1]-1]*y_mesh[(int)index_mesh[3*n+2]-1]-y_mesh[(int)index_mesh[3*n+1]-1]*x_mesh[(int)index_mesh[3*n+2]-1]
+		  +  x_mesh[(int)index_mesh[3*n+0]-1]*y_mesh[(int)index_mesh[3*n+1]-1]-y_mesh[(int)index_mesh[3*n+0]-1]*x_mesh[(int)index_mesh[3*n+1]-1]
+		  +  x_mesh[(int)index_mesh[3*n+2]-1]*y_mesh[(int)index_mesh[3*n+0]-1]-y_mesh[(int)index_mesh[3*n+2]-1]*x_mesh[(int)index_mesh[3*n+0]-1];
+
+
+		/*Go through x_grid and y_grid and interpolate if necessary*/
+		for (i=0;i<nlines;i++){
+
+			//exit if y not between y_tria_min and y_tria_max
+			if((y_grid[i]>y_tria_max) || (y_grid[i]<y_tria_min)) continue;
+
+			for(j=0;j<ncols; j++){
+
+				//exit if x not between x_tria_min and x_tria_max
+				if((x_grid[j]>x_tria_max) || (x_grid[j]<x_tria_min)) continue;
+
+				/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
+				area_1=((x_grid[j]-x_mesh[(int)index_mesh[3*n+2]-1])*(y_mesh[(int)index_mesh[3*n+1]-1]-y_mesh[(int)index_mesh[3*n+2]-1]) 
+							-  (y_grid[i]-y_mesh[(int)index_mesh[3*n+2]-1])*(x_mesh[(int)index_mesh[3*n+1]-1]-x_mesh[(int)index_mesh[3*n+2]-1]))/area;
+				/*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
+				area_2=((x_mesh[(int)index_mesh[3*n+0]-1]-x_mesh[(int)index_mesh[3*n+2]-1])*(y_grid[i]-y_mesh[(int)index_mesh[3*n+2]-1]) 
+							- (y_mesh[(int)index_mesh[3*n+0]-1]-y_mesh[(int)index_mesh[3*n+2]-1])*(x_grid[j]-x_mesh[(int)index_mesh[3*n+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_mesh[(int)index_mesh[3*n+0]-1]+area_2*data_mesh[(int)index_mesh[3*n+1]-1]+area_3*data_mesh[(int)index_mesh[3*n+2]-1];
+					}
+					else{
+						/*element interpolation*/
+						data_value=data_mesh[n];
+					}
+					if (isnan(data_value)) data_value=default_value;
+
+					/*insert value and go to the next point*/
+					griddata[i*ncols+j]=data_value;
+				}
+			}
+		}
+	}
+	if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
+
+	/*Assign output pointers:*/
+	*pgriddata=griddata;
+	*px_m=x_m;
+	*py_m=y_m;
+}
Index: /issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.h
===================================================================
--- /issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.h	(revision 2284)
+++ /issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.h	(revision 2284)
@@ -0,0 +1,12 @@
+/*!\file InterpFromMeshToGridx.h
+ * \brief: header file for Data interpolation routines.
+ */
+
+#ifndef _INTERPFROMMESHTOGRIDX_H
+#define _INTERPFROMMESHTOGRIDX_H
+
+#include "../toolkits/toolkits.h"
+
+int InterpFromMeshToGridx(double** px_m,double** py_m,double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh, int data_length, double cornereast,double cornernorth,double xposting,double yposting,int nlines,int ncols,double default_value);
+
+#endif /* _INTERPFROMMESHTOGRIDX_H */