Index: /issm/trunk-jpl/src/jl/md/utils.jl
===================================================================
--- /issm/trunk-jpl/src/jl/md/utils.jl	(revision 26711)
+++ /issm/trunk-jpl/src/jl/md/utils.jl	(revision 26712)
@@ -53,5 +53,63 @@
 	return output
 end# }}}
-function InterpFromMeshToMesh2d(index::Array,x::Vector,y::Vector,data::Vector,xout::Vector,yout::Vector) #{{{
+function InterpFromMeshToMesh2d(index_data::Array,x_data::Vector,y_data::Vector,data::Vector,xout::Vector,yout::Vector,default::Float64=NaN) #{{{
+
+	#Allocate output
+	nods_out = length(xout)
+	data_out = default*ones(nods_out)
+
+	#Interpolation type
+	data_length = size(data,1)
+	nods_data   = length(x_data)
+	nels_data   = size(index_data,1)
+	if(data_length==nods_data)
+		interpolation_type=1;
+	elseif (data_length==nels_data)
+		interpolation_type=2
+	else
+		error("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!")
+	end
+	xmin = minimum(xout); xmax = maximum(xout)
+	ymin = minimum(yout); ymax = maximum(yout)
+
+	for i in 1:nels_data
+
+		#skip element if no overlap
+		if (minimum(x_data[index_data[i,:]]) > xmax) continue end
+		if (minimum(y_data[index_data[i,:]]) > ymax) continue end
+		if (maximum(x_data[index_data[i,:]]) < xmin) continue end
+		if (maximum(y_data[index_data[i,:]]) < ymin) continue end
+
+		#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[index_data[i,2]]*y_data[index_data[i,3]]-y_data[index_data[i,2]]*x_data[index_data[i,3]] 
+				  +  x_data[index_data[i,1]]*y_data[index_data[i,2]]-y_data[index_data[i,1]]*x_data[index_data[i,2]] 
+				  +  x_data[index_data[i,3]]*y_data[index_data[i,1]]-y_data[index_data[i,3]]*x_data[index_data[i,1]])
+
+		for j in 1:nods_out
+			#Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area
+			area_1=((xout[j]-x_data[index_data[i,3]])*(y_data[index_data[i,2]]-y_data[index_data[i,3]])
+					 -  (yout[j]-y_data[index_data[i,3]])*(x_data[index_data[i,2]]-x_data[index_data[i,3]]))/area
+			#Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area
+			area_2=((x_data[index_data[i,1]]-x_data[index_data[i,3]])*(yout[j]-y_data[index_data[i,3]])
+					  - (y_data[index_data[i,1]]-y_data[index_data[i,3]])*(xout[j]-x_data[index_data[i,3]]))/area
+			#Get third area coordinate = 1-area1-area2
+			area_3=1-area_1-area_2
+
+			if (area_1>=0 && area_2>=0 && area_3>=0)
+				if (interpolation_type==1)
+					#nodal interpolation
+					data_out[j]=area_1*data[index_data[i,1]]+area_2*data[index_data[i,2]]+area_3*data[index_data[i,3]];
+				else
+					#element interpolation
+					data_out[j]=data[i];
+				end
+			end
+		end
+	end
+	return data_out
+
+	#OLD STUFF!!! not working...
+	
 	#prepare input arrays
 	nods = Cint(length(x))
@@ -126,2 +184,7 @@
 	end
 end #}}}
+function meshgrid(x::Vector, y::Vector)
+    X = [i for i in x, j in 1:length(y)]
+    Y = [j for i in 1:length(x), j in y]
+    return X, Y
+end