Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 12975)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 12976)
@@ -4305,5 +4305,5 @@
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  dp[NDOF2];
-	IssmDouble  vx,vy;
+	IssmDouble  vx,vy,vel;
 	GaussTria  *gauss                    = NULL;
 
@@ -4334,4 +4334,7 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
+		vel = sqrt(vx*vx+vy*vy);
+		vx  = vx/(vel+1.e-9);
+		vy  = vy/(vel+1.e-9);
 
 		/*J = 1/2 ( vx*dH/dx + vy*dH/dy )^2 */
@@ -4354,5 +4357,5 @@
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  dp[NDOF2];
-	IssmDouble  vx,vy;
+	IssmDouble  vx,vy,vel;
 	GaussTria  *gauss                    = NULL;
 
@@ -4383,4 +4386,7 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
+		vel = sqrt(vx*vx+vy*vy);
+		vx  = vx/(vel+1.e-9);
+		vy  = vy/(vel+1.e-9);
 
 		/*J = 1/2 ( -vy*dH/dx + vx*dH/dy )^2 */
@@ -4453,5 +4459,5 @@
 	IssmDouble  dbasis[NDOF2][NUMVERTICES];
 	IssmDouble  dH[2];
-	IssmDouble  v[2];
+	IssmDouble  vx,vy,vel;
 	GaussTria *gauss     = NULL;
 	int       *responses = NULL;
@@ -4498,13 +4504,19 @@
 			case ThicknessAlongGradientEnum:
 				weights_input->GetInputValue(&weight, gauss,resp);
-				vx_input->GetInputValue(&v[0],gauss);
-				vy_input->GetInputValue(&v[1],gauss);
-				for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*v[0]+dH[1]*v[1])*(dbasis[0][i]*v[0]+dbasis[1][i]*v[1])*Jdet*gauss->weight;
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vel = sqrt(vx*vx+vy*vy);
+				vx  = vx/(vel+1.e-9);
+				vy  = vy/(vel+1.e-9);
+				for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0][i]*vx+dbasis[1][i]*vy)*Jdet*gauss->weight;
 				break;
 			case ThicknessAcrossGradientEnum:
 				weights_input->GetInputValue(&weight, gauss,resp);
-				vx_input->GetInputValue(&v[0],gauss);
-				vy_input->GetInputValue(&v[1],gauss);
-				for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*(-v[1])+dH[1]*v[0])*(dbasis[0][i]*(-v[1])+dbasis[1][i]*v[0])*Jdet*gauss->weight;
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vel = sqrt(vx*vx+vy*vy);
+				vx  = vx/(vel+1.e-9);
+				vy  = vy/(vel+1.e-9);
+				for(i=0;i<numdof;i++) pe->values[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0][i]*(-vy)+dbasis[1][i]*vx)*Jdet*gauss->weight;
 				break;
 			default: