Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25935)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25936)
@@ -4140,8 +4140,7 @@
 		mk = 1.0 /3.0;
 		element->ElementSizes(&hx,&hy,&hz);
-		hk = max(max(hx, hy), hz);
-		// if(!element->IsOnDirichletBoundary()) {
-			Tau = - mk * hk * hk * 0.125 / viscosity;
-
+		// hk = max(max(hx, hy), hz);
+		hk = max(hx, hy);
+		Tau = - mk * hk * hk * 0.125 / viscosity;
 
 		for (int q=0;q<dim;q++){
@@ -4149,5 +4148,5 @@
 			for (int p=0;p<dim;p++){
 				for(int i=0;i<vnumnodes;i++){
-					SU[q*vnumnodes*(dim+1)+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
+					SU[q*numdof+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
 				}
 			}
@@ -4155,10 +4154,10 @@
 			for(int i=0;i<vnumnodes;i++){
 				for (int p=0;p<dim;p++){
-					SU[q*vnumnodes*(dim+1)+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
+					SU[q*numdof+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
 				}
 			}
 			/* column 4 for the pressure */
-			for(int i=0;i<vnumnodes;i++){
-					SU[q*vnumnodes*(dim+1)+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
+			for(int i=0;i<pnumnodes;i++){
+					SU[q*numdof+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
 			}
 		}
@@ -4226,4 +4225,5 @@
 	int vnumnodes = element->NumberofNodesVelocity();
 	int pnumnodes = element->NumberofNodesPressure();
+	int numdof    = vnumnodes*dim + pnumnodes;
 
 	/*Prepare coordinate system list*/
@@ -4291,11 +4291,7 @@
 		mk = 1.0 /3.0;
 		element->ElementSizes(&hx,&hy,&hz);
-		hk = max(max(hx, hy), hz);
-		// if(!element->IsOnDirichletBoundary()) {
-			Tau = - mk * hk * hk * 0.125 / viscosity;
-		// }
-		// else {
-		// 	Tau = 0.0;
-		// }
+		// hk = max(max(hx, hy), hz);
+		hk = max(hx, hy);
+		Tau = - mk * hk * hk * 0.125 / viscosity;
 
 		for (int q=0;q<dim;q++){
@@ -4303,5 +4299,5 @@
 			for (int p=0;p<dim;p++){
 				for(int i=0;i<vnumnodes;i++){
-					SU[q*vnumnodes*(dim+1)+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
+					SU[q*numdof+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
 				}
 			}
@@ -4309,13 +4305,12 @@
 			for(int i=0;i<vnumnodes;i++){
 				for (int p=0;p<dim;p++){
-					SU[q*vnumnodes*(dim+1)+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
+					SU[q*numdof+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
 				}
 			}
 			/* column 4 for the pressure */
-			for(int i=0;i<vnumnodes;i++){
-					SU[q*vnumnodes*(dim+1)+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
+			for(int i=0;i<pnumnodes;i++){
+					SU[q*numdof+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
 			}
 		}
-
 		for(i=0;i<vnumnodes;i++){
 			pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
@@ -4326,6 +4321,6 @@
 		}
 
-		for(int i=0;i<vnumnodes*(dim+1);i++){
-			pe->values[i] += Tau*Jdet*gauss->weight*rho_ice*(-gravity)*SU[(dim-1)*(dim+1)*vnumnodes+i];
+		for(int i=0;i<numdof;i++){
+			pe->values[i] += Tau*Jdet*gauss->weight*rho_ice*(-gravity)*SU[(dim-1)*numdof+i];
 		}
 	}
@@ -5159,5 +5154,5 @@
 		for(int i=0;i<vnumnodes;i++){
 			for(int d=0;d<dim;d++){
-				nsigma[d*vnumnodes+i] = 2.0 * viscosity * bed_normal[d]*vdbasisn[i];
+				nsigma[d*vnumnodes+i] = 1.0 * viscosity * bed_normal[d]*vdbasisn[i];
 			}
 		}
@@ -5172,5 +5167,5 @@
 					for (int q=0; q<dim; q++) {
 						/* gamma*(n*u)*(n*v) */
-						Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (-gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
+						Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
 						/* sigma(u)*(n*v) */
 						Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* nsigma[p*vnumnodes+i] * bed_normal[q] * vbasis[j];