Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 20658)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 20659)
@@ -482,4 +482,8 @@
 	int domaintype,dim=2;
 
+	/*Get approximation*/
+	int approximation;
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+
 	/* Get node coordinates and dof list: */
 	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
@@ -487,5 +491,4 @@
 	/*Retrieve all inputs we will be needing: */
 	this->FindParam(&domaintype,DomainTypeEnum);
-	if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
 	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
@@ -498,5 +501,17 @@
 		/*Compute strain rate and viscosity: */
 		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
-		this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
+		switch(approximation){
+			case SSAApproximationEnum:
+				this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
+				break;
+			case HOApproximationEnum:
+				this->material->ViscosityHO(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
+				break;
+			case FSApproximationEnum:
+				this->material->ViscosityFS(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input,NULL);
+				break;
+			default:
+				_error_("not supported yet");
+		}
 
 		/*Compute Stress*/