Index: /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp	(revision 20650)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp	(revision 20651)
@@ -385,4 +385,45 @@
 	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
 	ec_input->GetInputValue(&Ec,gauss);
+	es_input->GetInputValue(&Es,gauss);
+	ko_input->GetInputValue(&ko,gauss);
+
+	/*Compute viscosity*/
+	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+}/*}}}*/
+void  Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void  Matbafl::ResetHooks(){/*{{{*/
+
+	this->element=NULL;
+
+	/*Get Element type*/
+	this->helement->reset();
+
+}
+/*}}}*/
+IssmDouble Matbafl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble vorticity[3],vorticity_norm;
+	IssmDouble nrsp[3],nrsp_norm;
+	IssmDouble eps[3][3],epso;
+	IssmDouble epsprime[3],epsprime_norm;
+	IssmDouble E,lambdas;
+
+	/*Create vorticity vector*/
+	_assert_(dvx && dvy && dvz);
+	vorticity[0] =  dvz[1] - dvy[2];
+	vorticity[1] =  dvx[2] - dvz[0];
+	vorticity[2] =  dvy[0] - dvx[1];
+
+	/*Normalize*/
+	vorticity_norm = sqrt(vorticity[0]*vorticity[0] + vorticity[1]*vorticity[1] + vorticity[2]*vorticity[2]);
+	if(vorticity_norm==0){
+		vorticity[0] = 0.;
+		vorticity[1] = 0.;
+		vorticity[2] = 1.;
+	}
 	else{
 		vorticity[0] =vorticity[0]/vorticity_norm;