Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18832)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18833)
@@ -438,4 +438,58 @@
 	 }
 
+	/*Get viscosity*/
+	material->GetViscosity_B(&dmudB,eps_eff);
+
+	/*Assign output pointer*/
+	*pdmudB=dmudB;
+
+}
+/*}}}*/
+void       Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble dmudB;
+	IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
+	IssmDouble epsilon2d[2];/* epsilon=[exx,eyy,exy];    */
+	IssmDouble eps_eff;
+	IssmDouble eps0=1.e-27;
+
+	if(dim==3){
+		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
+		this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
+	}
+	else{
+		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
+		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1] + eps0*eps0);
+	}
+	/*Get viscosity*/
+	material->GetViscosity_B(&dmudB,eps_eff);
+
+	/*Assign output pointer*/
+	*pdmudB=dmudB;
+
+}
+/*}}}*/
+void       Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble dmudB;
+	IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
+	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
+	IssmDouble eps_eff;
+	IssmDouble eps0=1.e-27;
+
+	if(dim==3){
+		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
+		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
+	}
+	else{
+		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
+		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
+	}
 	/*Get viscosity*/
 	material->GetViscosity_B(&dmudB,eps_eff);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18832)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 18833)
@@ -66,4 +66,6 @@
 		void       DeleteMaterials(void);
 		void       dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
+		void       dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
+		void       dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
 		void       dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		IssmDouble Divergence(void);