Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17092)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17093)
@@ -413,7 +413,8 @@
 	 */
 
-	IssmDouble epsilonvx[6];
-	IssmDouble epsilonvy[6];
-	IssmDouble epsilonvz[6];
+	/*Intermediaries*/
+	IssmDouble dvx[3];
+	IssmDouble dvy[3];
+	IssmDouble dvz[3];
 
 	/*Check that both inputs have been found*/
@@ -423,10 +424,14 @@
 
 	/*Get strain rate assuming that epsilon has been allocated*/
-	vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
-	vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
-
-	/*Sum all contributions*/
-	for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
+	epsilon[0] = dvx[0];
+	epsilon[1] = dvy[1];
+	epsilon[2] = dvz[2];
+	epsilon[3] = 0.5*(dvx[1] + dvy[0]);
+	epsilon[4] = 0.5*(dvx[2] + dvz[0]);
+	epsilon[5] = 0.5*(dvy[2] + dvz[1]);
+
 }/*}}}*/
 void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
@@ -441,7 +446,7 @@
 	 */
 
-	int i;
-	IssmDouble epsilonvx[5];
-	IssmDouble epsilonvy[5];
+	/*Intermediaries*/
+	IssmDouble dvx[3];
+	IssmDouble dvy[3];
 
 	/*Check that both inputs have been found*/
@@ -451,15 +456,18 @@
 
 	/*Get strain rate assuming that epsilon has been allocated*/
-	vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
-
-	/*Sum all contributions*/
-	for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	epsilon[0] = dvx[0];
+	epsilon[1] = dvy[1];
+	epsilon[2] = 0.5*(dvx[1] + dvy[0]);
+	epsilon[3] = 0.5*dvx[2];
+	epsilon[4] = 0.5*dvy[2];
+
 }/*}}}*/
 void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
-	int i;
-	IssmDouble epsilonvx[3];
-	IssmDouble epsilonvy[3];
+	/*Intermediaries*/
+	IssmDouble dvx[3];
+	IssmDouble dvy[3];
 
 	/*Check that both inputs have been found*/
@@ -469,9 +477,10 @@
 
 	/*Get strain rate assuming that epsilon has been allocated*/
-	vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
-
-	/*Sum all contributions*/
-	for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+	epsilon[0] = dvx[0];
+	epsilon[1] = dvy[1];
+	epsilon[2] = 0.5*(dvx[1] + dvy[0]);
+
 }/*}}}*/
 void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/