Index: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4837)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4838)
@@ -3024,5 +3024,4 @@
 	int          numberofdofspernode;
 
-
 	/* 3d gaussian points: */
 	int     num_gauss,ig;
@@ -3097,4 +3096,10 @@
 	bool shelf;
 
+	/*inputs: */
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vxold_input=NULL;
+	Input* vyold_input=NULL;
+
 	/*retrieve inputs :*/
 	inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
@@ -3112,5 +3117,4 @@
 	  bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
 	  the stiffness matrix. */
-
 
 	if ((collapse==1) && (onbed==0)){
@@ -3136,5 +3140,4 @@
 		GetDofList(&doflist[0],&numberofdofspernode);
 
-
 		/*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
 		  get tria gaussian points as well as segment gaussian points. For tria gaussian 
@@ -3148,4 +3151,10 @@
 		GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
 
+		/*Retrieve all inputs we will be needing: */
+		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: */
 		for (ig1=0; ig1<num_area_gauss; ig1++){
@@ -3165,6 +3174,6 @@
 
 				/*Get strain rate from velocity: */
-				inputs->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum);
-				inputs->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum);
+				this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
+				this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
 
 				/*Get viscosity: */
@@ -3344,4 +3353,9 @@
 	bool isstokes;
 
+	/*inputs: */
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vz_input=NULL;
+
 	/*retrive parameters: */
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
@@ -3379,4 +3393,9 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
+
+	/*Retrieve all inputs we will be needing: */
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+	vz_input=inputs->GetInput(VzEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -3393,5 +3412,5 @@
 
 			/*Compute strain rate: */
-			inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
+			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
 
 			/*Get viscosity: */
@@ -3462,5 +3481,5 @@
 
 			/*Compute strain rate: */
-			inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
+			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
 
 			/*Get viscosity at last iteration: */
@@ -4094,4 +4113,6 @@
 	bool collapse;
 	bool onbed;
+	Input* surface_input=NULL;
+	Input* thickness_input=NULL;
 
 	/*retrieve inputs :*/
@@ -4135,4 +4156,8 @@
 		GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
 
+		/*Retrieve all inputs we will be needing: */
+		thickness_input=inputs->GetInput(ThicknessEnum);
+		surface_input=inputs->GetInput(SurfaceEnum);
+
 		/* Start  looping on the number of gaussian points: */
 		for (ig1=0; ig1<num_area_gauss; ig1++){
@@ -4150,8 +4175,8 @@
 
 				/*Compute thickness at gaussian point: */
-				inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum);
+				thickness_input->GetParameterValue(&thickness, gauss_coord);
 
 				/*Compute slope at gaussian point: */
-				inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum);
+				surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord);
 
 				/* Get Jacobian determinant: */
@@ -4293,4 +4318,8 @@
 	bool shelf;
 	bool isstokes;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vz_input=NULL;
+	Input* bed_input=NULL;
 
 	/*retrieve inputs :*/
@@ -4326,4 +4355,10 @@
 	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
+
+	/*Retrieve all inputs we will be needing: */
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+	vz_input=inputs->GetInput(VzEnum);
+	bed_input=inputs->GetInput(BedEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -4340,5 +4375,5 @@
 
 			/*Compute strain rate and viscosity: */
-			inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
+			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
 			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
@@ -4421,5 +4456,5 @@
 
 			/* Get bed at gaussian point */
-			inputs->GetParameterValue(&bed, gauss_coord,BedEnum);
+			bed_input->GetParameterValue(&bed, gauss_coord);
 
 			/*Get L matrix : */
@@ -4514,4 +4549,6 @@
 	bool onwater;
 	bool onbed;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
 
 	/*retrieve inputs :*/
@@ -4541,4 +4578,8 @@
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
 
+	/*Retrieve all inputs we will be needing: */
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+
 	/* Start  looping on the number of gaussian points: */
 	for (ig1=0; ig1<num_area_gauss; ig1++){
@@ -4557,6 +4598,6 @@
 			/*Get velocity derivative, with respect to x and y: */
 
-			inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum);
-			inputs->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord,VyEnum);
+			vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord);
+			vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord);
 			dudx=du[0];
 			dvdy=dv[1];
@@ -4725,4 +4766,7 @@
 	bool onbed;
 	bool shelf;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vz_input=NULL;
 
 	/*retrieve inputs :*/
@@ -4757,4 +4801,9 @@
 	/*Get gaussian points: */
 	GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
+
+	/*Retrieve all inputs we will be needing: */
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+	vz_input=inputs->GetInput(VzEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -4771,5 +4820,5 @@
 
 			/*Compute strain rate and viscosity: */
-			inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
+			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
 			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
@@ -6014,4 +6063,63 @@
 }
 /*}}}*/
+/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/
+void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){
+	/*Compute the 3d Blatter/PattynStrain Rate (5 components):
+	 *
+	 * epsilon=[exx eyy exy exz eyz]
+	 *
+	 * with exz=1/2 du/dz
+	 *      eyz=1/2 dv/dz
+	 *
+	 * the contribution of vz is neglected
+	 */
+
+	int i;
+
+	double epsilonvx[5];
+	double epsilonvy[5];
+
+	/*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->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
+	vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
+
+	/*Sum all contributions*/
+	for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+
+}
+/*}}}*/
+/*FUNCTION Penta::GetStrainRate3d{{{1*/
+void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
+	/*Compute the 3d Strain Rate (6 components):
+	 *
+	 * epsilon=[exx eyy ezz exy exz eyz]
+	 */
+
+	int i;
+
+	double epsilonvx[6];
+	double epsilonvy[6];
+	double epsilonvz[6];
+
+	/*Check that both inputs have been found*/
+	if (!vx_input || !vy_input || !vz_input){
+		ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input);
+	}
+
+	/*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(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
+
+}
+/*}}}*/
 /*FUNCTION Penta::GradjB {{{1*/
 void  Penta::GradjB(Vec gradient){
Index: issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.h	(revision 4837)
+++ issm/trunk/src/c/objects/Elements/Penta.h	(revision 4838)
@@ -168,6 +168,6 @@
 		void	  GetSolutionFromInputsDiagnosticVert(Vec solutiong);
 		void	  GetSolutionFromInputsThermal(Vec solutiong);
-		void	  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
-		void	  GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
+		void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
+		void    GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
 		Penta*  GetUpperElement(void);
 		void	  InputExtrude(int enum_type,bool only_if_collapsed);
