Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16421)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16422)
@@ -10693,6 +10693,5 @@
 
 	/*Initialize values*/
-	IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
-	IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
+	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
 	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
 	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
@@ -10701,4 +10700,9 @@
 	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);
@@ -10706,15 +10710,15 @@
 
 	/*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]];
+	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(&vvalues[0],nodes,vnumnodes,XYZEnum);
+	TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
 	for(i=0;i<vnumnodes;i++){
-		vx[i] = vvalues[i*NDOF3+0];
-		vy[i] = vvalues[i*NDOF3+1];
-		vz[i] = vvalues[i*NDOF3+2];
+		vx[i] = values[i*NDOF3+0];
+		vy[i] = values[i*NDOF3+1];
+		vz[i] = values[i*NDOF3+2];
 		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
@@ -10722,5 +10726,5 @@
 	}
 	for(i=0;i<pnumnodes;i++){
-		pressure[i] = pvalues[i];
+		pressure[i] = values[vnumdof+i];
 		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
 	}
@@ -10728,6 +10732,6 @@
 	/*Recondition pressure and compute vel: */
 	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
-	for(i=0;i<vnumnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i = 0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
+	for(i = 0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
 	/*Now, we have to move the previous inputs  to old 
@@ -10751,8 +10755,8 @@
 	xDelete<IssmDouble>(vy);
 	xDelete<IssmDouble>(vx);
-	xDelete<IssmDouble>(vvalues);
-	xDelete<IssmDouble>(pvalues);
+	xDelete<IssmDouble>(values);
 	xDelete<int>(vdoflist);
 	xDelete<int>(pdoflist);
+	xDelete<int>(cs_list);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16421)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16422)
@@ -4417,4 +4417,5 @@
 	GetDofListPressure(&pdoflist,GsetEnum);
 
+	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
 	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];