Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 22522)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 22523)
@@ -156,9 +156,9 @@
 	/*Intermediaries */
 	int  stabilization,dim, domaintype, calvinglaw;
-	int i, row, col;
+	int i,j,k,row, col;
 	IssmDouble kappa;
 	IssmDouble Jdet, dt, D_scalar;
 	IssmDouble h,hx,hy,hz;
-	IssmDouble vel;
+	IssmDouble vel,v[3],w[3],c[3],m[3],dlsf[3];
 	IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
 	IssmDouble calvingmax, calvinghaf, heaviside, haf_eps;
@@ -184,12 +184,9 @@
 	ElementMatrix* Ke       = basalelement->NewElementMatrix();
 	IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
-	IssmDouble*    v        = xNew<IssmDouble>(dim);
-	IssmDouble*    w        = xNew<IssmDouble>(dim);
-	IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
-	IssmDouble*    m        = xNewZeroInit<IssmDouble>(dim);
-	IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
+	IssmDouble*    dbasis   = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    Bprime = NULL;
+	if(stabilization==2){
+		Bprime   = xNew<IssmDouble>(dim*numnodes);
+	}
 
 	/*Retrieve all inputs and parameters*/
@@ -283,19 +280,19 @@
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		D_scalar=gauss->weight*Jdet;
 
 		/* Transient */
 		if(dt!=0.){
-			basalelement->NodalFunctions(basis,gauss);
-			TripleMultiply(basis,numnodes,1,0,
-						&D_scalar,1,1,0,
-						basis,1,numnodes,0,
-						&Ke->values[0],1);
-			D_scalar*=dt;
+			for(i=0;i<numnodes;i++){
+				for(j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
+				}
+			}
+			D_scalar=D_scalar*dt;
 		}
 
 		/* Advection */
-		GetB(B,basalelement,xyz_list,gauss); 
-		GetBprime(Bprime,basalelement,xyz_list,gauss); 
 		vx_input->GetInputValue(&v[0],gauss);
 		vy_input->GetInputValue(&v[1],gauss); 
@@ -458,17 +455,11 @@
 
 		/*Compute D*/
-		for(row=0;row<dim;row++){
-			for(col=0;col<dim;col++){
-				if(row==col)
-				 D[row*dim+col]=D_scalar*w[row];
-				else
-				 D[row*dim+col]=0.;
+		for(i=0;i<numnodes;i++){
+			for(j=0;j<numnodes;j++){
+				for(k=0;k<dim;k++){
+					Ke->values[i*numnodes+j] += D_scalar*w[k]*dbasis[k*numnodes+j]*basis[i];
+				}
 			}
 		}
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					Bprime,dim,numnodes,0,
-					&Ke->values[0],1);
 
 		/* Stabilization */
@@ -485,15 +476,11 @@
 				h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) ); 
 				kappa=h*vel/2.;
-				for(row=0;row<dim;row++)
-					for(col=0;col<dim;col++)
-					if(row==col)
-						D[row*dim+col]=D_scalar*kappa;
-					else
-						D[row*dim+col]=0.;
-
-				TripleMultiply(Bprime,dim,numnodes,1,
-							D,dim,dim,0,
-							Bprime,dim,numnodes,0,
-							&Ke->values[0],1);
+				for(i=0;i<numnodes;i++){
+					for(j=0;j<numnodes;j++){
+						for(k=0;k<dim;k++){
+							Ke->values[i*numnodes+j] += D_scalar*kappa*dbasis[k*numnodes+j]*dbasis[k*numnodes+i];
+						}
+					}
+				}
 				break;	
 			case 2:
@@ -501,10 +488,12 @@
 				basalelement->ElementSizes(&hx,&hy,&hz);
 				h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
+				IssmDouble D[3*3];
 				for(row=0;row<dim;row++) 
 					for(col=0;col<dim;col++) 
 						D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
+				GetBprime(Bprime,basalelement,xyz_list,gauss);
 
 				TripleMultiply(Bprime,dim,numnodes,1,
-							D,dim,dim,0,
+							&D[0],dim,dim,0,
 							Bprime,dim,numnodes,0,
 							&Ke->values[0],1);
@@ -518,12 +507,6 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(dbasis);
 	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(v);
-	xDelete<IssmDouble>(w);
-	xDelete<IssmDouble>(c);
-	xDelete<IssmDouble>(m);
-	xDelete<IssmDouble>(dlsf);
 	delete gauss;
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 22522)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 22523)
@@ -461,5 +461,4 @@
 		}
 		else if(stabilization==2){
-			element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 			tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
 			for(int i=0;i<numnodes;i++){