Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18699)
@@ -1106,4 +1106,6 @@
 				name==EplHeadOldEnum ||
 				name==StressIntensityFactorEnum ||
+				name==StrainRateparallelEnum ||
+				name==StrainRateperpendicularEnum ||
 				name==HydrologydcEplThicknessOldEnum ||
 				name==HydrologydcEplInitialThicknessEnum ||
@@ -1208,4 +1210,12 @@
 				input=this->inputs->GetInput(output_enum);
 				break;
+			case StrainRateparallelEnum:
+				this->StrainRateparallel();
+				input=this->inputs->GetInput(output_enum);
+				break;
+			case StrainRateperpendicularEnum:
+				this->StrainRateperpendicular();
+				input=this->inputs->GetInput(output_enum);
+				break;
 			default:
 				_error_("input "<<EnumToStringx(output_enum)<<" not found in element");
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18699)
@@ -220,4 +220,6 @@
 		virtual void   ComputeDeviatoricStressTensor(void)=0;
 		virtual void	StressIntensityFactor(void)=0;
+		virtual void	StrainRateparallel(void)=0;
+		virtual void   StrainRateperpendicular(void)=0;
 
 		virtual void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18699)
@@ -422,4 +422,93 @@
 	this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
 	this->InputExtrude(StressIntensityFactorEnum,-1);
+}
+/*}}}*/
+void       Penta::StrainRateparallel(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainparallel[NUMVERTICES];
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	
+	/*Retrieve all inputs we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
+	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
+	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainxx_input->GetInputValue(&strainxx,gauss);
+		strainxy_input->GetInputValue(&strainxy,gauss);
+		strainyy_input->GetInputValue(&strainyy,gauss);
+
+		/*strainparallel= Strain rate along the ice flow direction */
+		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
+void       Penta::StrainRateperpendicular(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainperpendicular[NUMVERTICES];
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
+	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
+	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
+
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainxx_input->GetInputValue(&strainxx,gauss);
+		strainxy_input->GetInputValue(&strainxy,gauss);
+		strainyy_input->GetInputValue(&strainyy,gauss);
+
+		/*strainperpendicular= Strain rate perpendicular to the ice flow direction */
+		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18699)
@@ -58,4 +58,6 @@
 		void   ComputeDeviatoricStressTensor();
 		void   StressIntensityFactor();
+		void   StrainRateparallel();
+		void   StrainRateperpendicular();
 		void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18699)
@@ -58,4 +58,6 @@
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
 		void        StressIntensityFactor(void){_error_("not implemented yet");};
+		void			StrainRateparallel(void){_error_("not implemented yet");};
+		void			StrainRateperpendicular(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18699)
@@ -58,4 +58,6 @@
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
 		void        StressIntensityFactor(void){_error_("not implemented yet");};
+		void        StrainRateparallel(void){_error_("not implemented yet");};
+		void        StrainRateperpendicular(void){_error_("not implemented yet");};
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18699)
@@ -271,4 +271,93 @@
 	this->inputs->AddInput(new TriaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum));
 	this->inputs->AddInput(new TriaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
+void       Tria::StrainRateperpendicular(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainperpendicular[NUMVERTICES];
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
+	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
+	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
+
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainxx_input->GetInputValue(&strainxx,gauss);
+		strainxy_input->GetInputValue(&strainxy,gauss);
+		strainyy_input->GetInputValue(&strainyy,gauss);
+
+		/*strainperpendicular= Strain rate perpendicular to the ice flow direction */
+		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
+void       Tria::StrainRateparallel(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainparallel[NUMVERTICES];
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	
+	/*Retrieve all inputs we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input* strainxx_input=inputs->GetInput(StrainRatexxEnum);             _assert_(strainxx_input);
+	Input* strainxy_input=inputs->GetInput(StrainRatexyEnum);             _assert_(strainxy_input);
+	Input* strainyy_input=inputs->GetInput(StrainRateyyEnum);             _assert_(strainyy_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainxx_input->GetInputValue(&strainxx,gauss);
+		strainxy_input->GetInputValue(&strainxy,gauss);
+		strainyy_input->GetInputValue(&strainyy,gauss);
+
+		/*strainparallel= Strain rate along the ice flow direction */
+		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
 
 	/*Clean up and return*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18698)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18699)
@@ -57,4 +57,6 @@
 		void        ComputeSurfaceNormalVelocity();
 		void        StressIntensityFactor(void){_error_("not implemented yet");};
+		void			StrainRateparallel();
+		void			StrainRateperpendicular();
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
 		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);