Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25610)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 25611)
@@ -2728,11 +2728,8 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = basalelement->GetNumberOfNodes();
-	int numdof   = numnodes*2;
 
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = basalelement->NewElementMatrix(MLHOApproximationEnum);
-	IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
+	IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -2749,15 +2746,24 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBSSA(B,basalelement,2,xyz_list,gauss_base);
-		this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
 
 		element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
 
-		for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet;
-
-		TripleMultiply(B,3,numdof,1,
-					D,3,3,0,
-					Bprime,3,numdof,0,
-					&Ke->values[0],1);
+		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*viscosity*(
+							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*viscosity*(
+							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*viscosity*(
+							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*viscosity*(
+							dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+							);
+			}
+		}
 	}
 
@@ -2770,7 +2776,5 @@
 	basalelement->DeleteMaterials(); delete basalelement;
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(dbasis);
 	return Ke;
 }/*}}}*/