Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6220)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6221)
@@ -2962,6 +2962,8 @@
 		D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
 
-		MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
-		MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0);
+		TripleMultiply(&B_conduct[0][0],3,numdof,1,
+					&D[0][0],3,3,0,
+					&B_conduct[0][0],3,numdof,0,
+					&Ke_gaussian_conduct[0][0],0);
 
 		/*Advection: */
@@ -2981,6 +2983,8 @@
 		D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
 
-		MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
-		MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
+		TripleMultiply(&B_advec[0][0],3,numdof,1,
+					&D[0][0],3,3,0,
+					&Bprime_advec[0][0],3,numdof,0,
+					&Ke_gaussian_advec[0][0],0);
 
 		/*Transient: */
@@ -2991,6 +2995,8 @@
 			D_scalar=D_scalar;
 
-			MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
-			MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0);
+			TripleMultiply(&L[0],numdof,1,0,
+						&D_scalar,1,1,0,
+						&L[0],1,numdof,0,
+						&Ke_gaussian_transient[0][0],0);
 		}
 		else{
@@ -3009,6 +3015,8 @@
 			GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
 
-			MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
-			MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
+			TripleMultiply(&B_artdiff[0][0],2,numdof,1,
+						&K[0][0],2,2,0,
+						&B_artdiff[0][0],2,numdof,0,
+						&Ke_gaussian_artdiff[0][0],0);
 		}
 		else if(artdiff==2){
@@ -3575,6 +3583,8 @@
 		for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
 
-		MatrixMultiply(&B[0][0],8,numdofbubble,1,&D[0][0],8,8,0,&tBD[0][0],0);
-		MatrixMultiply(&tBD[0][0],numdofbubble,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0);
+		TripleMultiply(&B[0][0],8,numdofbubble,1,
+					&D[0][0],8,8,0,
+					&B_prime_bubble[0][0],8,3,0,
+					&Ke_gaussian[0][0],0);
 
 		for(i=0;i<numdofbubble;i++) for(j=0;j<NDOF3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
@@ -5578,6 +5588,8 @@
 
 	/*Multiply matrices to create the reduce matrix Ke_reduced */
-	MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0);
-	MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Kbi[0][0],3,24,0,&Kright[0][0],0);
+	TripleMultiply(&Kib[0][0],24,3,0,
+				&Kbbinv[0][0],3,3,0,
+				&Kbi[0][0],3,24,0,
+				&Kright[0][0],0);
 
 	/*Affect value to the reduced matrix */
@@ -5622,6 +5634,7 @@
 
 	/*Multiply matrices to create the reduce matrix Ke_reduced */
-	MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0);
-	MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Pb[0],3,1,0,&Pright[0],0);
+	TripleMultiply(&Kib[0][0],24,3,0,
+				&Kbbinv[0][0],3,3,0,
+				&Pb[0],3,1,0,&Pright[0],0);
 
 	/*Affect value to the reduced matrix */
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 6220)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 6221)
@@ -3164,6 +3164,8 @@
 		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
 
-		MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
-		MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
+		TripleMultiply(&L[0],numdof,1,0,
+					&D_scalar,1,1,0,
+					&L[0],1,numdof,0,
+					&Ke_gaussian[0][0],0);
 
 		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
@@ -3482,6 +3484,8 @@
 		if(dt) D_scalar=dt*D_scalar;
 
-		MatrixMultiply(&l1l2l3[0],numdof,1,0,&D_scalar,1,1,0,&tl1l2l3D[0],0);
-		MatrixMultiply(&tl1l2l3D[0],numdof,1,0,&l1l2l3[0],1,numdof,0,&Ke_gaussian[0][0],0);
+		TripleMultiply(&l1l2l3[0],numdof,1,0,
+					&D_scalar,1,1,0,
+					&l1l2l3[0],1,numdof,0,
+					&Ke_gaussian[0][0],0);
 
 		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];