Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26131)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 26132)
@@ -2737,5 +2737,5 @@
 
 	/*Intermediaries*/
-	IssmDouble  viscosity,Jdet;
+	IssmDouble  viscosity,Jdet,thickness,n;
 	IssmDouble *xyz_list = NULL;
 
@@ -2748,11 +2748,14 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = basalelement->NewElementMatrix(MLHOApproximationEnum);
-	IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes); // like SSA
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes); // like SSA
 
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* vx_input      = element->GetInput(VxEnum);      _assert_(vx_input);
-	Input* vy_input      = element->GetInput(VyEnum);      _assert_(vy_input);
+	Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* surface_input  = element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input);
+	Input* n_input			 = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -2764,21 +2767,69 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
-
+		element->NodalFunctions(basis, gauss);
+
+		thickness_input->GetInputValue(&thickness, gauss);
+		n_input->GetInputValue(&n,gauss);
 		element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
-
-		for(int i=0;i<numnodes;i++){
+		
+		for(int i=0;i<numnodes;i++){//shape functions on tria element
 			for(int j=0;j<numnodes;j++){
-				Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
+				Ke->values[(i+0)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
 							4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
-							);
-				Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
+							);//K11
+				Ke->values[(i+0)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+							)*(n+1)/(n+2);//K12
+				Ke->values[(i+0)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
 							2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
-							);
-				Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
+							);//K13
+				Ke->values[(i+0)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+                     2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
+                     )*(n+1)/(n+2);//K14
+				
+				Ke->values[(i+3)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+							)*(n+1)/(n+2);//K21
+				Ke->values[(i+3)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+							)*2*pow(n+1,2)/((2*n+3)*(n+2)) 
+							+ 
+							gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
+							;//K22
+				Ke->values[(i+3)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+							2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
+							)*(n+1)/(n+2);//K23
+				Ke->values[(i+3)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+							2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
+							)*2*pow(n+1,2)/((2*n+3)*(n+2));//K24
+				
+				Ke->values[(i+6)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
 							2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
-							);
-				Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
-							dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
-							);
+							);//K31
+				Ke->values[(i+6)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+							2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
+							)*(n+1)/(n+2);//K32
+				Ke->values[(i+6)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
+							);//K33
+				Ke->values[(i+6)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
+							)*(n+1)/(n+2);//K34
+
+				Ke->values[(i+9)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                     2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
+                     )*(n+1)/(n+2);//K41
+				Ke->values[(i+9)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                     2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
+                     )*2*pow(n+1,2)/((2*n+3)*(n+2));//K42
+				Ke->values[(i+9)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
+							)*(n+1)/(n+2);//K43
+				Ke->values[(i+9)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+							4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
+							)*2*pow(n+1,2)/((2*n+3)*(n+2))
+							+ 
+							gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
+							;//K44
 			}
 		}
@@ -2794,4 +2845,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(basis);
 	return Ke;
 }/*}}}*/