Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 25264)
+++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 25265)
@@ -62,5 +62,6 @@
 	/*Initialize Element vector and other vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -73,25 +74,24 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		GetB(B,element,dim,xyz_list,gauss);
-		D=gauss->weight*Jdet*dt;
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		/*vertical diffusion*/
-		TripleMultiply(B,1,numnodes,1,
-					&D,1,1,0,
-					B,1,numnodes,0,
-					&Ke->values[0],1);
+		D=gauss->weight*Jdet*dt;
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D*dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j];
+			}
+		}
 
-		/*Next value*/
+		/*Next value (Transient)*/
 		D=gauss->weight*Jdet;
-		element->NodalFunctions(B,gauss);
-		TripleMultiply(B,numnodes,1,0,
-					&D,1,1,0,
-					B,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[j]*basis[i];
 	} 
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	return Ke;
@@ -137,24 +137,4 @@
 	return pe;
 }/*}}}*/
-void           DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
-		where hi is the interpolation function for node i.*/
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[i] = dbasis[(dim-1)*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
 void           DepthAverageAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
Index: /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h	(revision 25264)
+++ /issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.h	(revision 25265)
@@ -26,5 +26,4 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void           GetB(IssmDouble* B,Element* element,int dim,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/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 25264)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 25265)
@@ -446,9 +446,5 @@
 		if(dt!=0.){
 			D_scalar=gauss->weight*Jdet;
-			for(int i=0;i<numnodes;i++){
-				for(int j=0;j<numnodes;j++){
-					Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
-				}
-			}
+			for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
 			D_scalar=D_scalar*dt;
 		}