Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18240)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18241)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 24
+//#define FSANALYTICAL 12
 
 /*Model processing*/
@@ -3246,6 +3246,6 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-	IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+	IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
+	IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -3271,9 +3271,9 @@
 		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
-
-		TripleMultiply(B,dim-1,numdof,1,
-					D,dim-1,dim-1,0,
-					B,dim-1,numdof,0,
+		for(int i=0;i<dim;i++) D[i*dim+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
+
+		TripleMultiply(B,dim,numdof,1,
+					D,dim,dim,0,
+					B,dim,numdof,0,
 					&Ke->values[0],1);
 	}
@@ -3394,6 +3394,6 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-	IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+	IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
+	IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
 
 	/*Retrieve all inputs and parameters*/
@@ -3431,9 +3431,9 @@
 		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
-
-		TripleMultiply(B,dim-1,numdof,1,
-					D,dim-1,dim-1,0,
-					B,dim-1,numdof,0,
+		for(int i=0;i<dim;i++) D[i*dim+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
+
+		TripleMultiply(B,dim,numdof,1,
+					D,dim,dim,0,
+					B,dim,numdof,0,
 					&Ke->values[0],1);
 	}
@@ -4429,4 +4429,8 @@
 			B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
 			B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
+
+			B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
+			B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
+			B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
 		}
 		for(int i=0;i<pnumnodes;i++){
@@ -4437,6 +4441,9 @@
 	else{
 		for(int i=0;i<vnumnodes;i++){
-			B[2*i+0] = vbasis[i];
-			B[2*i+1] = 0.;
+			B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
+			B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
+
+			B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
+			B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
 		}