Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4752)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4753)
@@ -2940,4 +2940,9 @@
 
 	/*inputs: */
+	Input* thickness_input=NULL;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vxold_input=NULL;
+	Input* vyold_input=NULL;
 	bool onwater,shelf;
 
@@ -2958,4 +2963,11 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/*Retrieve all inputs we will be needing: */
+	thickness_input=inputs->GetInput(ThicknessEnum);
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+	vxold_input=inputs->GetInput(VxOldEnum);
+	vyold_input=inputs->GetInput(VyOldEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2969,9 +2981,9 @@
 
 		/*Compute thickness at gaussian point: */
-		inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum);
+		thickness_input->GetParameterValue(&thickness, gauss_l1l2l3);
 
 		/*Get strain rate from velocity: */
-		inputs->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum);
-		inputs->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxOldEnum,VyOldEnum);
+		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,vx_input,vy_input);
+		this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,vxold_input,vyold_input);
 
 		/*Get viscosity: */
@@ -3073,4 +3085,5 @@
 	bool shelf;
 	int  drag_type;
+	Input* surface_input=NULL;
 
 	/*retrive parameters: */
@@ -3080,4 +3093,5 @@
 	inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
 	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	surface_input=inputs->GetInput(SurfaceEnum);
 	
 	/* Get node coordinates and dof list: */
@@ -3114,8 +3128,7 @@
 		GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
 
-
 		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
 		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
-		inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum);
+		surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
 		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
 
@@ -4212,8 +4225,14 @@
 	bool onwater;
 	int  drag_type;
+	Input* thickness_input=NULL;
+	Input* surface_input=NULL;
+	Input* drag_input=NULL;
 
 	/*retrieve inputs :*/
 	inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
 	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	thickness_input=inputs->GetInput(ThicknessEnum);
+	surface_input=inputs->GetInput(SurfaceEnum);
+	drag_input=inputs->GetInput(DragCoefficientEnum);
 
 	/*First, if we are on water, return empty vector: */
@@ -4237,11 +4256,11 @@
 
 		/*Compute thickness at gaussian point: */
-		inputs->GetParameterValue(&thickness, &gauss_l1l2l3[0],ThicknessEnum);
-		inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum);
+		thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]);
+		surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
 		
 		/*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 
 		 * element itself: */
 		if(drag_type==1){
-			inputs->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0],DragCoefficientEnum);
+			drag_input->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0]);
 		}
 
@@ -5473,4 +5492,30 @@
 	/*Add value to global vector*/
 	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
+
+}
+/*}}}*/
+/*FUNCTION Tria::GetStrainRate2d{{{1*/
+void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){
+	/*Compute the 2d Strain Rate (3 components):
+	 *
+	 * epsilon=[exx eyy exy]
+	 */
+
+	int i;
+
+	double epsilonvx[3];
+	double epsilonvy[3];
+
+	/*Check that both inputs have been found*/
+	if (!vx_input || !vy_input){
+		ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
+	}
+
+	/*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];
 
 }
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4752)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4753)
@@ -153,9 +153,10 @@
 		void	  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
 		void	  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
-		void    GetParameterValue(double* pvalue,Node* node,int enumtype);
-		void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
+		void      GetParameterValue(double* pvalue,Node* node,int enumtype);
+		void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
 		void	  GetSolutionFromInputsDiagnosticHoriz(Vec solution);
 		void	  GetSolutionFromInputsAdjointHoriz(Vec solution);
 		void	  GetSolutionFromInputsDiagnosticHutter(Vec solution);
+		void      GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
 		void	  GradjDragStokes(Vec gradient);
 		void	  InputUpdateFromSolutionAdjointHoriz( double* solution);