Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16372)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16373)
@@ -8186,6 +8186,4 @@
 	int         analysis_type,approximation;
 	IssmDouble  alpha2,Jdet2d;
-	IssmDouble  FSreconditioning,viscosity;
-	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -8213,5 +8211,4 @@
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -8230,7 +8227,4 @@
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
 		GetLFS(BFriction,gauss);
-
-		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16372)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16373)
@@ -3208,9 +3208,7 @@
 	/*Intermediaries */
 	int         i,j;
-	int         analysis_type,approximation;
+	int         analysis_type;
 	int         indices[2];
 	IssmDouble  alpha2,Jdet;
-	IssmDouble  FSreconditioning,viscosity;
-	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_seg[NUMVERTICES1D][3];
@@ -3232,11 +3230,10 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke        = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
-	IssmDouble*    BFriction = xNew<IssmDouble>(2*numdof);
-	IssmDouble*    D         = xNewZeroInit<IssmDouble>(2*2);
+	IssmDouble*    BFriction = xNew<IssmDouble>(1*numdof);
+	IssmDouble     D;
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -3256,6 +3253,13 @@
 
 		GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
-
-		_error_("STOP");
+		GetLFS(BFriction,gauss);
+
+		friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
+		D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx
+
+		TripleMultiply(BFriction,1,numdof,1,
+					&D,1,1,0,
+					BFriction,1,numdof,0,
+					Ke->values,1);
 	}
 
@@ -3267,5 +3271,4 @@
 	delete friction;
 	xDelete<IssmDouble>(BFriction);
-	xDelete<IssmDouble>(D);
 	xDelete<int>(cs_list);
 	return Ke;
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16372)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16373)
@@ -469,4 +469,37 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetLFS{{{*/
+void TriaRef::GetLFS(IssmDouble* LFS, GaussTria* gauss){
+	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For node i, Li can be expressed in the actual coordinate system
+	 * by: 
+	 *       Li=[ h 0 0 ]
+	 * where h is the interpolation function for node i.
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int pnumnodes = this->NumberofNodesPressure();
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumdof   = pnumnodes;
+	int vnumdof   = vnumnodes*NDOF2;
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
+	GetNodalFunctionsVelocity(vbasis,gauss);
+
+	/*Build LFS: */
+	for(int i=0;i<vnumnodes;i++){
+		LFS[2*i+0] = vbasis[i];
+		LFS[2*i+1] = 0.;
+	}
+
+	for(int i=0;i<pnumnodes;i++){
+		LFS[i+vnumdof+0] = 0.;
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(vbasis);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16372)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16373)
@@ -37,4 +37,5 @@
 		void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
+		void GetLFS(IssmDouble* LFS, GaussTria* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int finiteelement);
Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 16372)
+++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 16373)
@@ -91,5 +91,9 @@
 	if (Neff<0)Neff=0;
 
-	if(strcmp(element_type,"2d")==0){
+	if(strcmp(element_type,"1d")==0){
+		this->GetInputValue(&vx, gauss,vxenum);
+		vmag=sqrt(vx*vx);
+	}
+	else if(strcmp(element_type,"2d")==0){
 		this->GetInputValue(&vx, gauss,vxenum);
 		this->GetInputValue(&vy, gauss,vyenum);