Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22512)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 22513)
@@ -1307,6 +1307,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum);
-	IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
-	IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -1341,12 +1340,22 @@
 		}
 
-		this->GetBSSAFriction(B,element,dim,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
 		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
-		for(int i=0;i<dim;i++) D[i*dim+i]=alpha2*gauss->weight*Jdet;
-
-		TripleMultiply(B,dim,numdof,1,
-					D,dim,dim,0,
-					B,dim,numdof,0,
-					&Ke->values[0],1);
+
+		if(dim==2){
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[2*i*2*numnodes+2*j]       += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+					Ke->values[(2*i+1)*2*numnodes+2*j+1] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+				}
+			}
+		}
+		else{
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+				}
+			}
+		}
 	}
 
@@ -1358,6 +1367,5 @@
 	delete friction;
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(basis);
 	return Ke;
 }/*}}}*/
@@ -1434,8 +1442,7 @@
 
 	/*Intermediaries*/
-	int         dim,domaintype,bsize;
+	int         dim,domaintype;
 	IssmDouble  viscosity,newviscosity,oldviscosity;
 	IssmDouble  viscosity_overshoot,thickness,Jdet;
-	IssmDouble  D_scalar;
 	IssmDouble *xyz_list = NULL;
 
@@ -1443,7 +1450,7 @@
 	element->FindParam(&domaintype,DomainTypeEnum);
 	switch(domaintype){
-		case Domain2DverticalEnum:   dim = 1; bsize = 1; break;
-		case Domain2DhorizontalEnum: dim = 2; bsize = 3; break;
-		case Domain3DEnum:           dim = 2; bsize = 3; break;
+		case Domain2DverticalEnum:   dim = 1; break;
+		case Domain2DhorizontalEnum: dim = 2; break;
+		case Domain3DEnum:           dim = 2; break;
 		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 	}
@@ -1455,7 +1462,5 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
-	IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -1478,19 +1483,38 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBSSA(B,element,dim,xyz_list,gauss);
-		this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss);
-
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		thickness_input->GetInputValue(&thickness, gauss);
 		element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
 		element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
-		thickness_input->GetInputValue(&thickness, gauss);
-
 		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
-		for(int i=0;i<bsize;i++) D[i*bsize+i]=D_scalar;
-
-		TripleMultiply(B,bsize,numdof,1,
-					D,bsize,bsize,0,
-					Bprime,bsize,numdof,0,
-					&Ke->values[0],1);
+
+		if(dim==2){
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*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*newviscosity*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*newviscosity*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*newviscosity*thickness*(
+								dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+								);
+				}
+			}
+		}
+		else{
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += gauss->weight*Jdet*newviscosity*thickness*(
+								4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
+								);
+				}
+			}
+		}
 	}
 
@@ -1501,7 +1525,5 @@
 	delete gauss;
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(dbasis);
 	return Ke;
 }/*}}}*/