Index: /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp	(revision 24935)
+++ /issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp	(revision 24936)
@@ -164,6 +164,6 @@
 	/*Initialize Element vector and other vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*		dbasis = xNew<IssmDouble>(2*numnodes);
 	IssmDouble     D[2][2];
 
@@ -190,6 +190,6 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		GetB(B,element,xyz_list,gauss);
-		GetBprime(Bprime,element,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -201,19 +201,12 @@
 		D_scalar=gauss->weight*Jdet;
 
-		D[0][0]=D_scalar*dvxdx;
-		D[0][1]=0.;
-		D[1][0]=0.;
-		D[1][1]=D_scalar*dvydy;
-		TripleMultiply(B,2,numnodes,1,
-					&D[0][0],2,2,0,
-					B,2,numnodes,0,
-					&Ke->values[0],1);
-
-		D[0][0]=D_scalar*vx;
-		D[1][1]=D_scalar*vy;
-		TripleMultiply(B,2,numnodes,1,
-					&D[0][0],2,2,0,
-					Bprime,2,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				/*\phi_i \phi_j \nabla\cdot v*/
+				Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+				/*\phi_i v\cdot\nabla\phi_j*/
+				Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
+			}
+		}
 
 		if(stabilization==1){
@@ -239,8 +232,13 @@
 			D[0][1]=D_scalar*D[0][1];
 			D[1][1]=D_scalar*D[1][1];
-			TripleMultiply(Bprime,2,numnodes,1,
-						&D[0][0],2,2,0,
-						Bprime,2,numnodes,0,
-						&Ke->values[0],1);
+
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += (
+								dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) +
+								dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j]) 
+								);
+				}
+			}
 		}
 	}
@@ -248,6 +246,6 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	return Ke;
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 24935)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 24936)
@@ -216,8 +216,5 @@
 		D = -gauss->weight*Jdet*normal[2];
 
-		TripleMultiply( basis,1,numnodes,1,
-					&D,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
 	}
 
@@ -256,8 +253,5 @@
 		D = -gauss->weight*Jdet*normal[2];
 
-		TripleMultiply( basis,1,numnodes,1,
-					&D,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D*basis[i]*basis[j];
 	}