Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 25262)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 25263)
@@ -454,6 +454,5 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 	IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
 
@@ -461,5 +460,4 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
-	//printf("dt %f\n", dt);
 	element->FindParam(&stabilization,DamageStabilizationEnum);
 	Input2* vx_input = element->GetInput2(VxEnum); _assert_(vx_input);
@@ -479,4 +477,5 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
@@ -490,35 +489,33 @@
 		}
 
+		/*Transient term*/
 		D_scalar=gauss->weight*Jdet;
-		TripleMultiply(basis,1,numnodes,1,
-					&D_scalar,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
-
-		GetB(B,element,dim,xyz_list,gauss);
-		GetBprime(Bprime,element,dim,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];
 
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
-		if(dim==3) dvzdz=dvz[2];
 		D_scalar=dt*gauss->weight*Jdet;
-
-		D[0*dim+0]=D_scalar*dvxdx;
-		D[1*dim+1]=D_scalar*dvydy;
-		if(dim==3) D[2*dim+2]=D_scalar*dvzdz;
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					B,dim,numnodes,0,
-					&Ke->values[0],1);
-
-		D[0*dim+0]=D_scalar*vx;
-		D[1*dim+1]=D_scalar*vy;
-		if(dim==3) D[2*dim+2]=D_scalar*vz;
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					Bprime,dim,numnodes,0,
-					&Ke->values[0],1);
+		if(dim==2){
+			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]);
+				}
+			}
+		}
+		else{/*3D*/
+			_assert_(dim==3);
+			dvzdz=dvz[2];
+			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+dvzdz);
+					/*\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] + vz*dbasis[2*numnodes+j]);
+				}
+			}
+		}
 
 		if(stabilization==2){
@@ -566,4 +563,12 @@
 				D[0*dim+1]=D_scalar*D[0*dim+1];
 				D[1*dim+1]=D_scalar*D[1*dim+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*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
+									dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
+									);
+					}
+				}
 			}
 			else if(dim==3){
@@ -577,11 +582,15 @@
 				D[1*dim+2]=D_scalar*D[1*dim+2];
 				D[2*dim+2]=D_scalar*D[2*dim+2];
-			}
-			TripleMultiply(Bprime,dim,numnodes,1,
-						D,dim,dim,0,
-						Bprime,dim,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*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j] + D[0*dim+2]*dbasis[2*numnodes+j]) +
+									dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j] + D[1*dim+2]*dbasis[2*numnodes+j]) +
+									dbasis[2*numnodes+i] *(D[2*dim+0]*dbasis[0*numnodes+j] + D[2*dim+1]*dbasis[1*numnodes+j] + D[2*dim+2]*dbasis[2*numnodes+j]) 
+									);
+					}
+				}
+			}
+		}
 	}
 
@@ -589,6 +598,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(D);
 	delete gauss;
@@ -666,61 +674,4 @@
 	return pe;
 }/*}}}*/
-void           DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,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++){
-		for(int j=0;j<dim;j++){
-			B[numnodes*j+i] = basis[i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,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>(dim*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		for(int j=0;j<dim;j++){
-			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/
 void           DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum);
@@ -788,7 +739,6 @@
 	/*Initialize Element vector and other vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -803,24 +753,22 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		GetB(B,element,dim,xyz_list,gauss);
-		GetBprime(Bprime,element,dim,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
 		vxaverage_input->GetInputValue(&vx,gauss);
 		vyaverage_input->GetInputValue(&vy,gauss);
 
-		D[0*dim+0] = -gauss->weight*vx*Jdet;
-		D[1*dim+1] = -gauss->weight*vy*Jdet;
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					Bprime,dim,numnodes,0,
-					&Ke->values[0],1);
-
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				/*\phi_i v\cdot\nabla\phi_j*/
+				Ke->values[i*numnodes+j] += -gauss->weight*Jdet*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
+			}
+		}
 	}
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
 	delete gauss;
 	return Ke;
@@ -854,8 +802,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];
 	}
 
Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 25262)
+++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h	(revision 25263)
@@ -30,6 +30,4 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBprime(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/L2ProjectionBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp	(revision 25262)
+++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp	(revision 25263)
@@ -108,9 +108,5 @@
 		basalelement->NodalFunctions(basis,gauss);
 		D=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/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 25262)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 25263)
@@ -337,9 +337,5 @@
 		/*Transient term*/
 		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[i]*basis[j];
-			}
-		}
+		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*/