Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18183)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18184)
@@ -2855,7 +2855,6 @@
 	if(fe_FS==XTaylorHoodEnum)
 	 Ke1=CreateKMatrixFSViscousXTH(element);
-	else if(fe_FS==LATaylorHoodEnum){
+	else if(fe_FS==LATaylorHoodEnum)
 	 Ke1=CreateKMatrixFSViscousLATH(element);
-	}
 	else
 	 Ke1=CreateKMatrixFSViscous(element);
@@ -2902,5 +2901,5 @@
 	IssmDouble*    B        = xNew<IssmDouble>(epssize*numdof);
 	IssmDouble*    Bprime   = xNew<IssmDouble>(epssize*numdof);
-	IssmDouble*    BtBUzawa = xNew<IssmDouble>(numdof*pnumdof);
+	IssmDouble*    BtBUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
 	IssmDouble*    BU       = xNew<IssmDouble>(pnumdof);
 	IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
@@ -2920,6 +2919,6 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBFS(B,element,dim,xyz_list,gauss);
-		this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
+		this->GetBFSvel(B,element,dim,xyz_list,gauss);
+		this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss);
 
 		element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
@@ -4121,4 +4120,131 @@
 	xDelete<IssmDouble>(pbasis);
 }/*}}}*/
+void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 
+	 * For node i, Bvi can be expressed in the actual coordinate system
+	 * by: 	   Bvi=[ dphi/dx          0        ]
+	 *					 [   0           dphi/dy     ]
+	 *					 [ 1/2*dphi/dy    1/2*dphi/dx]
+	 *
+	 *
+	 *	In 3d:
+	 *     	   Bvi=[ dh/dx          0             0      ]
+	 *					[   0           dh/dy           0      ]
+	 *					[   0             0           dh/dz    ]
+	 *					[ 1/2*dh/dy    1/2*dh/dx        0      ]
+	 *					[ 1/2*dh/dz       0         1/2*dh/dx  ]
+	 *					[   0          1/2*dh/dz    1/2*dh/dy  ]
+	 *
+	 *	where phi is the finiteelement function for node i.
+	 *	Same thing for Bb except the last column that does not exist.
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
+	element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+
+	/*Build B: */
+	if(dim==2){
+		for(int i=0;i<vnumnodes;i++){
+			B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*0+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+		}
+	}
+	else{
+		for(int i=0;i<vnumnodes;i++){
+			B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*0+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*0+dim*i+2] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*1+dim*i+2] = 0.;
+			B[(dim*vnumnodes)*2+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*2+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*3+dim*i+2] = 0.;
+			B[(dim*vnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*4+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*5+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
+		}
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(vdbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
+	 *	For node i, Bi' can be expressed in the actual coordinate system
+	 *	by: 
+	 *			Bvi' = [  dphi/dx     0     ]
+	 *					 [     0      dphi/dy ]
+	 *					 [  dphi/dy   dphi/dx ]
+	 *
+	 *	In 3d
+	 *     	   Bvi=[ dh/dx     0        0    ]
+	 *					[   0      dh/dy      0    ]
+	 *					[   0        0      dh/dz  ]
+	 *					[ dh/dy    dh/dx      0    ]
+	 *					[ dh/dz      0      dh/dx  ]
+	 *					[   0      dh/dz    dh/dy  ]
+	 *	where phi is the finiteelement function for node i.
+	 *	In 3d:
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
+	element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+
+	/*Build B_prime: */
+	if(dim==2){
+		for(int i=0;i<vnumnodes;i++){
+			Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
+		}
+	}
+	else{
+		for(int i=0;i<vnumnodes;i++){
+			Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes)*0+dim*i+2] = 0.;
+			Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*1+dim*i+2] = 0.;
+			Bprime[(dim*vnumnodes)*2+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes)*2+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*3+dim*i+2] = 0.;
+			Bprime[(dim*vnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*4+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*5+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
+		}
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(vdbasis);
+}/*}}}*/
 void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. 
@@ -4284,5 +4410,5 @@
 	/*Add value to global vector*/
 	solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
-	solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
+	if(pnumdof>0) solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
 
 	/*Free ressources:*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18183)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18184)
@@ -81,4 +81,6 @@
 		void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);