Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 22510)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 22511)
@@ -569,7 +569,5 @@
 	IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis   = xNew<IssmDouble>(3*numnodes);
-	IssmDouble*    B        = xNew<IssmDouble>(3*numnodes);
 	IssmDouble*    Bprime   = xNew<IssmDouble>(3*numnodes);
-	IssmDouble     D[3][3]  = {0.};
 	IssmDouble     K[3][3];
 
@@ -600,39 +598,39 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
 		D_scalar=gauss->weight*Jdet;
 		if(dt!=0.) D_scalar=D_scalar*dt;
 
 		/*Conduction: */
-		GetBConduct(B,element,xyz_list,gauss); 
-		D[0][0]=D_scalar*kappa/rho_ice;
-		D[1][1]=D_scalar*kappa/rho_ice;
-		D[2][2]=D_scalar*kappa/rho_ice;
-		TripleMultiply(B,3,numnodes,1,
-					&D[0][0],3,3,0,
-					B,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] += D_scalar*kappa/rho_ice*(
+							dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+							);
+			}
+		}
 
 		/*Advection: */
-		GetBAdvec(B,element,xyz_list,gauss); 
-		GetBAdvecprime(Bprime,element,xyz_list,gauss); 
 		vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
 		vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
 		vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
-		D[0][0]=D_scalar*vx;
-		D[1][1]=D_scalar*vy;
-		D[2][2]=D_scalar*vz;
-		TripleMultiply(B,3,numnodes,1,
-					&D[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] += D_scalar*(
+							vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
+							);
+			}
+		}
 
 		/*Transient: */
 		if(dt!=0.){
 			D_scalar=gauss->weight*Jdet;
-			element->NodalFunctions(basis,gauss);
-			TripleMultiply(basis,numnodes,1,0,
-						&D_scalar,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_scalar*basis[j]*basis[i];
+				}
+			}
 			D_scalar=D_scalar*dt;
 		}
@@ -678,5 +676,4 @@
 	xDelete<IssmDouble>(basis);
 	xDelete<IssmDouble>(dbasis);
-	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
 	delete gauss;
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 22510)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 22511)
@@ -242,5 +242,5 @@
 
 	/*Intermediaries*/
-	IssmDouble  D,Jdet;
+	IssmDouble  Jdet;
 	IssmDouble *xyz_list = NULL;
 
@@ -250,6 +250,6 @@
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix(NoneApproximationEnum);
-	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -262,12 +262,14 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetB(B,element,xyz_list,gauss);
-		this->GetBprime(Bprime,element,xyz_list,gauss);
-		D=gauss->weight*Jdet;
-
-		TripleMultiply(B,1,numnodes,1,
-					&D,1,1,0,
-					Bprime,1,numnodes,0,
-					&Ke->values[0],1);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+							basis[j]*dbasis[2*numnodes+i]
+							);
+			}
+		}
 	}
 
@@ -275,6 +277,6 @@
 	delete gauss;
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(basis);
 	return Ke;
 
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 22510)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 22511)
@@ -368,7 +368,6 @@
 	int         stabilization;
 	IssmDouble  Jdet,dt,u,v,w,um,vm,wm,vel;
-	IssmDouble  h,hx,hy,hz,vx,vy,vz;
+	IssmDouble  h,hx,hy,hz,vx,vy,vz,D_scalar;
 	IssmDouble  tau_parameter,diameter;
-	IssmDouble  D_scalar;
 	IssmDouble* xyz_list = NULL;
 
@@ -380,7 +379,5 @@
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(3*numnodes);
 	IssmDouble*    Bprime = xNew<IssmDouble>(3*numnodes);
-	IssmDouble     D[3][3]={0.};
 	IssmDouble     K[3][3];
 
@@ -409,39 +406,39 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
 		D_scalar=gauss->weight*Jdet;
 		if(dt!=0.) D_scalar=D_scalar*dt;
 
 		/*Conduction: */
-		GetBConduct(B,element,xyz_list,gauss); 
-		D[0][0]=D_scalar*kappa;
-		D[1][1]=D_scalar*kappa;
-		D[2][2]=D_scalar*kappa;
-		TripleMultiply(B,3,numnodes,1,
-					&D[0][0],3,3,0,
-					B,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] += D_scalar*kappa*(
+							dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+							);
+			}
+		}
 
 		/*Advection: */
-		GetBAdvec(B,element,xyz_list,gauss); 
-		GetBAdvecprime(Bprime,element,xyz_list,gauss); 
 		vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
 		vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
 		vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
-		D[0][0]=D_scalar*vx;
-		D[1][1]=D_scalar*vy;
-		D[2][2]=D_scalar*vz;
-		TripleMultiply(B,3,numnodes,1,
-					&D[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] += D_scalar*(
+							vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
+							);
+			}
+		}
 
 		/*Transient: */
 		if(dt!=0.){
 			D_scalar=gauss->weight*Jdet;
-			element->NodalFunctions(basis,gauss);
-			TripleMultiply(basis,numnodes,1,0,
-						&D_scalar,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_scalar*basis[j]*basis[i];
+				}
+			}
 			D_scalar=D_scalar*dt;
 		}
@@ -456,5 +453,4 @@
 			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); 
 
@@ -486,5 +482,4 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
 	xDelete<IssmDouble>(basis);