Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18192)
@@ -113,4 +113,75 @@
 	/*Assign output pointer*/
 	*ptransform=transform;
+}
+/*}}}*/
+void       Element::ComputeStrainRate(){/*{{{*/
+
+	int         dim;
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[6];
+
+	/*Retrieve all inputs we will be needing: */
+	this->GetVerticesCoordinates(&xyz_list);
+	parameters->FindParam(&dim,DomainDimensionEnum);
+	Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input=NULL;
+	if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);}
+
+	/*Allocate arrays*/
+	int numvertices = this->GetNumberOfVertices();
+	IssmDouble* eps_xx = xNew<IssmDouble>(numvertices);
+	IssmDouble* eps_yy = xNew<IssmDouble>(numvertices);
+	IssmDouble* eps_zz = xNew<IssmDouble>(numvertices);
+	IssmDouble* eps_xy = xNew<IssmDouble>(numvertices);
+	IssmDouble* eps_xz = xNew<IssmDouble>(numvertices);
+	IssmDouble* eps_yz = xNew<IssmDouble>(numvertices);
+
+	/* Start looping on the number of vertices: */
+	Gauss* gauss=this->NewGauss();
+	for (int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Compute strain rate viscosity and pressure: */
+		if(dim==2)
+		 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		else
+		 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+
+		if(dim==2){
+			 /* epsilon=[exx,eyy,exy];*/
+			eps_xx[iv]=epsilon[0]; 
+			eps_yy[iv]=epsilon[1];
+			eps_xy[iv]=epsilon[2];
+		}
+		else{
+			/*epsilon=[exx eyy ezz exy exz eyz]*/
+			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->AddInput(StrainRatexxEnum,eps_xx,P1Enum);
+	this->AddInput(StrainRatexyEnum,eps_xy,P1Enum);
+	this->AddInput(StrainRatexzEnum,eps_xz,P1Enum);
+	this->AddInput(StrainRateyyEnum,eps_yy,P1Enum);
+	this->AddInput(StrainRateyzEnum,eps_yz,P1Enum);
+	this->AddInput(StrainRatezzEnum,eps_zz,P1Enum);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(eps_xx);
+	xDelete<IssmDouble>(eps_yy);
+	xDelete<IssmDouble>(eps_zz);
+	xDelete<IssmDouble>(eps_xy);
+	xDelete<IssmDouble>(eps_xz);
+	xDelete<IssmDouble>(eps_yz);
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18192)
@@ -60,4 +60,5 @@
 		/* bool       AnyActive(void); */
 		void       CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
+		void       ComputeStrainRate();
 		void       Echo();
 		void       DeepEcho();
@@ -209,6 +210,4 @@
 		virtual void   ComputeSigmaNN(void)=0;
 		virtual void   ComputeBasalStress(Vector<IssmDouble>* sigma_b)=0;
-		virtual void   ComputeStrainRate()=0;
-		virtual void   ComputeStrainRate(Vector<IssmDouble>* eps)=0;
 		virtual void   ComputeStressTensor(void)=0;
 		virtual void   ComputeDeviatoricStressTensor(void)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18192)
@@ -223,52 +223,4 @@
 	/*Add value to output*/
 	sigma_b->SetValue(id-1,value,INS_VAL);
-}
-/*}}}*/
-void       Penta::ComputeStrainRate(){/*{{{*/
-
-	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/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 18192)
@@ -54,6 +54,4 @@
 		IssmDouble CharacteristicLength(void){_error_("not implemented yet");};
 		void   ComputeBasalStress(Vector<IssmDouble>* sigma_b);
-		void   ComputeStrainRate();
-		void   ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
 		void   ComputeSigmaNN(){_error_("not implemented yet");};
 		void   ComputeStressTensor();
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 18192)
@@ -55,6 +55,4 @@
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
 		void        ComputeSigmaNN(){_error_("not implemented yet");};
-		void        ComputeStrainRate(){_error_("not implemented yet");};
-		void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
 		void        ComputeStressTensor(){_error_("not implemented yet");};
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 18192)
@@ -55,6 +55,4 @@
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b){_error_("not implemented yet");};
 		void        ComputeSigmaNN(){_error_("not implemented yet");};
-		void        ComputeStrainRate(){_error_("not implemented yet");};
-		void        ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
 		void        ComputeStressTensor(){_error_("not implemented yet");};
 		void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18192)
@@ -196,49 +196,4 @@
 		delete gauss;
 	}
-}
-/*}}}*/
-void       Tria::ComputeStrainRate(){/*{{{*/
-
-	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};
-	GaussTria*     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 GaussTria();
-	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;
-
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18191)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18192)
@@ -53,6 +53,4 @@
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
 		void        ComputeSigmaNN();
-		void		   ComputeStrainRate();
-		void		   ComputeStrainRate(Vector<IssmDouble>* eps){_error_("not implemented yet");};
 		void        ComputeStressTensor();
 		void        ComputeDeviatoricStressTensor();