Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 25265)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 25266)
@@ -161,6 +161,5 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*		dbasis = xNew<IssmDouble>(2*numnodes);
 	IssmDouble     D[2][2]={0.};
 
@@ -183,4 +182,5 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
@@ -189,31 +189,20 @@
 		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
 
+		/*Transient term*/
 		D_scalar=gauss->weight*Jdet;
-
-		TripleMultiply(basis,1,numnodes,1,
-					&D_scalar,1,1,0,
-					basis,1,numnodes,0,
-					Ke->values,1);
-
-		GetB(B,element,xyz_list,gauss);
-		GetBprime(Bprime,element,xyz_list,gauss);
-
+		for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
+
+		/*Advection terms*/
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
 		D_scalar=dt*gauss->weight*Jdet;
-
-		D[0][0]=D_scalar*dvxdx;
-		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]);
+			}
+		}
 
 		/*Artificial diffusivity*/
@@ -223,8 +212,12 @@
 		D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
 		D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
-		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]) 
+							);
+			}
+		}
 	}
 
@@ -232,6 +225,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	return Ke;
@@ -285,59 +277,4 @@
 	delete gauss;
 	return pe;
-}/*}}}*/
-void           HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ N ]
-	 *          [ N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	element->NodalFunctions(basis,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = basis[i];
-		B[numnodes*1+i] = basis[i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ dN/dx ]
-	 *                [ dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-
 }/*}}}*/
 void           HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h	(revision 25265)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.h	(revision 25266)
@@ -27,6 +27,4 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 25265)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 25266)
@@ -1013,8 +1013,5 @@
 
 		D=gauss->weight*Jdet;
-		TripleMultiply(basis,1,numnodes,1,
-					&D,1,1,0,
-					basis,1,numnodes,0,
-					&Me->values[0],1);
+		for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Me->values[i*numnodes+j] += D*basis[i]*basis[j];
 	}