Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16400)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16401)
@@ -2456,7 +2456,5 @@
 
 		gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
-
 		slope_input->GetInputValue(&slope,gauss);
-
 		IssmDouble theta = atan(slope);
 
@@ -4404,4 +4402,5 @@
 
 	/*Initialize values*/
+	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
 	IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
 	IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
@@ -4411,14 +4410,22 @@
 	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
 
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
 	/*Get dof list: */
 	GetDofListVelocity(&vdoflist,GsetEnum);
 	GetDofListPressure(&pdoflist,GsetEnum);
 
+	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
+	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
+
+	/*Transform solution in Cartesian Space*/
+	TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
+
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
-	for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
-
-	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYEnum);
+	for(i=0;i<vnumdof;i++) vvalues[i]=values[i];
+	for(i=0;i<pnumdof;i++) pvalues[i]=values[vnumnodes+i];
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
@@ -4456,8 +4463,10 @@
 	xDelete<IssmDouble>(vy);
 	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
 	xDelete<IssmDouble>(vvalues);
 	xDelete<IssmDouble>(pvalues);
 	xDelete<int>(vdoflist);
 	xDelete<int>(pdoflist);
+	xDelete<int>(cs_list);
 }
 /*}}}*/