Index: /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp	(revision 25260)
+++ /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp	(revision 25261)
@@ -109,8 +109,5 @@
 		D=latentheat/heatcapacity*gauss->weight*Jdet;
 
-		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];
 	}
 
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 25260)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 25261)
@@ -362,9 +362,5 @@
 		D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
 		if(reCast<bool,IssmDouble>(dt)) D=dt*D;
-		TripleMultiply(basis,numnodes,1,0,
-					&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];
 	}
 
@@ -395,5 +391,4 @@
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(3*numnodes);
 	IssmDouble     K[3][3];
 
@@ -468,10 +463,14 @@
 			K[2][0]=h/(2.*vel)*fabs(vz*vx);  K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
 			for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
-			GetBAdvecprime(Bprime,element,xyz_list,gauss);
-
-			TripleMultiply(Bprime,3,numnodes,1,
-						&K[0][0],3,3,0,
-						Bprime,3,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] *(K[0][0]*dbasis[0*numnodes+j] + K[0][1]*dbasis[1*numnodes+j]+ K[0][2]*dbasis[2*numnodes+j]) +
+								dbasis[1*numnodes+i] *(K[1][0]*dbasis[0*numnodes+j] + K[1][1]*dbasis[1*numnodes+j]+ K[1][2]*dbasis[2*numnodes+j]) +
+								dbasis[2*numnodes+i] *(K[2][0]*dbasis[0*numnodes+j] + K[2][1]*dbasis[1*numnodes+j]+ K[2][2]*dbasis[2*numnodes+j]) 
+								);
+				}
+			}
 		}
 		else if(stabilization==2){
@@ -513,5 +512,4 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(Bprime);
 	xDelete<IssmDouble>(basis);
 	xDelete<IssmDouble>(dbasis);
@@ -744,91 +742,4 @@
 	return pe;
 
-}/*}}}*/
-void           ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by:
-	 *       Bi_advec =[ h ]
-	 *                 [ h ]
-	 *                 [ h ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(1*NUMNODESP1)
-	 */
-
-	/*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];
-		B[numnodes*2+i] = basis[i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by:
-	 *       Biprime_advec=[ dh/dx ]
-	 *                     [ dh/dy ]
-	 *                     [ dh/dz ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = dbasis[0*numnodes+i];
-		B[numnodes*1+i] = dbasis[1*numnodes+i];
-		B[numnodes*2+i] = dbasis[2*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}/*}}}*/
-void           ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by:
-	 *       Bi_conduct=[ dh/dx ]
-	 *                  [ dh/dy ]
-	 *                  [ dh/dz ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = dbasis[0*numnodes+i];
-		B[numnodes*1+i] = dbasis[1*numnodes+i];
-		B[numnodes*2+i] = dbasis[2*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void           ThermalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h	(revision 25260)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h	(revision 25261)
@@ -31,7 +31,4 @@
 		ElementVector* CreatePVectorShelf(Element* element);
 		ElementVector* CreatePVectorVolume(Element* element);
-		void           GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBConduct(IssmDouble* B,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);