Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10277)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10278)
@@ -5289,17 +5289,4 @@
 	h=sqrt(2*this->GetArea());
 
-	///*Create Artificial diffusivity once for all if requested*/
-	//if(stabilization){
-	//	gauss=new GaussTria();
-	//	gauss->GaussCenter();
-	//	GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
-	//	delete gauss;
-
-	//	vxaverage_input->GetInputAverage(&vx);
-	//	vyaverage_input->GetInputAverage(&vy);
-	//	K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
-	//	K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
-	//}
-
 	/*Start looping on the number of gaussian points:*/
 	gauss=new GaussTria(2);
@@ -5338,4 +5325,5 @@
 
 		if(stabilization==1){
+			/*Streamline upwinding*/
 			vel=sqrt(pow(vx,2.)+pow(vy,2.));
 			K[0][0]=h/(2*vel)*vx*vx;
@@ -5343,12 +5331,19 @@
 			K[0][1]=h/(2*vel)*vx*vy;
 			K[1][1]=h/(2*vel)*vy*vy;
+		}
+		else if(stabilization==2){
+			/*MacAyeal*/
+			vxaverage_input->GetInputAverage(&vx);
+			vyaverage_input->GetInputAverage(&vy);
+			K[0][0]=h/2.0*fabs(vx) /sqrt(2*sqrt(3)); // the second part should not be there
+			K[0][1]=0.;
+			K[1][0]=0.;
+			K[1][1]=h/2.0*fabs(vy) /sqrt(2*sqrt(3)); // the second part should not be there
+		}
+		if(stabilization==1 || stabilization==2){
 			KDL[0][0]=DL_scalar*K[0][0];
 			KDL[1][0]=DL_scalar*K[1][0];
 			KDL[0][1]=DL_scalar*K[0][1];
 			KDL[1][1]=DL_scalar*K[1][1];
-
-			//KDL[0][0]=DL_scalar*K[0][0];
-			//KDL[1][1]=DL_scalar*K[1][1];
-
 			TripleMultiply( &Bprime[0][0],2,numdof,1,
 						&KDL[0][0],2,2,0,