Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17811)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17812)
@@ -279,6 +279,48 @@
 void  Penta::ComputeStrainRate(Vector<IssmDouble>* eps){
 
-	_error_("Not implemented yet");
-
+	IssmDouble      xyz_list[NUMVERTICES][3];
+	IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble      eps_xx[NUMVERTICES];
+	IssmDouble		 eps_yy[NUMVERTICES];
+	IssmDouble		 eps_zz[NUMVERTICES];
+	IssmDouble      eps_xy[NUMVERTICES];
+	IssmDouble		 eps_xz[NUMVERTICES];
+	IssmDouble		 eps_yz[NUMVERTICES];
+	GaussPenta*     gauss=NULL;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);             _assert_(vz_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Compute strain rate viscosity and pressure: */
+		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+
+		eps_xx[iv]=epsilon[0]; 
+		eps_yy[iv]=epsilon[1];
+		eps_zz[iv]=epsilon[2];
+		eps_xy[iv]=epsilon[3];
+		eps_xz[iv]=epsilon[4];
+		eps_yz[iv]=epsilon[5];
+	}
+
+	/*Add Stress tensor components into inputs*/
+	this->inputs->AddInput(new PentaInput(StrainRatexxEnum,&eps_xx[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StrainRatexyEnum,&eps_xy[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StrainRatexzEnum,&eps_xz[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StrainRateyyEnum,&eps_yy[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StrainRateyzEnum,&eps_yz[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StrainRatezzEnum,&eps_zz[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17811)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17812)
@@ -204,5 +204,46 @@
 /*FUNCTION Tria::ComputeStrainRate {{{*/
 void  Tria::ComputeStrainRate(Vector<IssmDouble>* eps){
-	_error_("Not Implemented yet");
+
+	IssmDouble      xyz_list[NUMVERTICES][3];
+	IssmDouble      epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble      eps_xx[NUMVERTICES];
+	IssmDouble		 eps_yy[NUMVERTICES];
+	IssmDouble		 eps_zz[NUMVERTICES]={0,0,0};
+	IssmDouble      eps_xy[NUMVERTICES];
+	IssmDouble		 eps_xz[NUMVERTICES]={0,0,0};
+	IssmDouble		 eps_yz[NUMVERTICES]={0,0,0};
+	GaussPenta*     gauss=NULL;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Compute strain rate viscosity and pressure: */
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+
+		eps_xx[iv]=epsilon[0]; 
+		eps_yy[iv]=epsilon[1];
+		eps_xy[iv]=epsilon[2];
+	}
+
+	/*Add Stress tensor components into inputs*/
+	this->inputs->AddInput(new TriaInput(StrainRatexxEnum,&eps_xx[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRatexyEnum,&eps_xy[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRatexzEnum,&eps_xz[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRateyyEnum,&eps_yy[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRateyzEnum,&eps_yz[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(StrainRatezzEnum,&eps_zz[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+
 }
 /*}}}*/
