Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22533)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22534)
@@ -3135,5 +3135,5 @@
 
 	/*Intermediaries*/
-	int         i,dim,epssize;
+	int         i,dim;
 	IssmDouble  viscosity,FSreconditioning,Jdet;
 	IssmDouble *xyz_list = NULL;
@@ -3141,6 +3141,4 @@
 	/*Get problem dimension*/
 	element->FindParam(&dim,DomainDimensionEnum);
-	if(dim==2) epssize = 3;
-	else       epssize = 6;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -3148,5 +3146,4 @@
 	int pnumnodes = element->NumberofNodesPressure();
 	int numdof    = vnumnodes*dim + pnumnodes;
-	int bsize     = epssize + 2;
 
 	/*Prepare coordinate system list*/
@@ -3157,8 +3154,7 @@
 
 	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
+	ElementMatrix* Ke   = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
+	IssmDouble* pbasis  = xNew<IssmDouble>(pnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -3167,5 +3163,5 @@
 	Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
 	Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-	Input* vz_input;
+	Input* vz_input = NULL;
 	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
 
@@ -3176,15 +3172,59 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBFS(B,element,dim,xyz_list,gauss);
-		this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
-
+		element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+		element->NodalFunctionsPressure(pbasis,gauss);
 		element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
-		for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
-		for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
-
-		TripleMultiply(B,bsize,numdof,1,
-					D,bsize,bsize,0,
-					Bprime,bsize,numdof,0,
-					&Ke->values[0],1);
+
+		if(dim==3){
+			/*Stress balance*/
+			for(int i=0;i<vnumnodes;i++){
+				for(int j=0;j<vnumnodes;j++){
+					Ke->values[(3*i+0)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
+					Ke->values[(3*i+0)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
+					Ke->values[(3*i+0)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
+					Ke->values[(3*i+1)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
+					Ke->values[(3*i+1)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
+					Ke->values[(3*i+1)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
+					Ke->values[(3*i+2)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
+					Ke->values[(3*i+2)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
+					Ke->values[(3*i+2)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + 2.*vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
+				}
+				for(int k=0;k<pnumnodes;k++){
+					Ke->values[(3*i+0)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
+					Ke->values[(3*i+1)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
+					Ke->values[(3*i+2)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[2*vnumnodes+i]);
+				}
+			}
+			/*Incompressibility*/
+			for(int k=0;k<pnumnodes;k++){
+				for(int j=0;j<vnumnodes;j++){
+					Ke->values[(3*vnumnodes+k)*numdof+3*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
+					Ke->values[(3*vnumnodes+k)*numdof+3*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
+					Ke->values[(3*vnumnodes+k)*numdof+3*j+2] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[2*vnumnodes+j]*pbasis[k]);
+				}
+			}
+		}
+		else{
+			/*Stress balance*/
+			for(int i=0;i<vnumnodes;i++){
+				for(int j=0;j<vnumnodes;j++){
+					Ke->values[(2*i+0)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
+					Ke->values[(2*i+0)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
+					Ke->values[(2*i+1)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
+					Ke->values[(2*i+1)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
+				}
+				for(int k=0;k<pnumnodes;k++){
+					Ke->values[(2*i+0)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
+					Ke->values[(2*i+1)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
+				}
+			}
+			/*Incompressibility*/
+			for(int k=0;k<pnumnodes;k++){
+				for(int j=0;j<vnumnodes;j++){
+					Ke->values[(2*vnumnodes+k)*numdof+2*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
+					Ke->values[(2*vnumnodes+k)*numdof+2*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
+				}
+			}
+		}
 	}
 
@@ -3195,7 +3235,6 @@
 	delete gauss;
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(pbasis);
+	xDelete<IssmDouble>(vdbasis);
 	xDelete<int>(cs_list);
 	return Ke;
@@ -3706,6 +3745,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
-	IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
+	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -3741,12 +3779,26 @@
 		}
 
-		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		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);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+
+		if(dim==3){
+			/*Stress balance*/
+			for(int i=0;i<vnumnodes;i++){
+				for(int j=0;j<vnumnodes;j++){
+					Ke->values[(3*i+0)*numdof+3*j+0] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
+					Ke->values[(3*i+1)*numdof+3*j+1] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
+					Ke->values[(3*i+2)*numdof+3*j+2] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
+				}
+			}
+		}
+		else{
+			/*Stress balance*/
+			for(int i=0;i<vnumnodes;i++){
+				for(int j=0;j<vnumnodes;j++){
+					Ke->values[(2*i+0)*numdof+2*j+0] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
+					Ke->values[(2*i+1)*numdof+2*j+1] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j];
+				}
+			}
+		}
 	}
 
@@ -3757,6 +3809,5 @@
 	delete friction;
 	xDelete<IssmDouble>(xyz_list_base);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(vbasis);
 	return Ke;
 }/*}}}*/