Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10276)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10277)
@@ -501,6 +501,6 @@
 	int        stabilization;
 	int        i,j,ig,dim;
-	double     Jdettria,DL_scalar,dt;
-	double     vx,vy,dvxdx,dvydy;
+	double     Jdettria,DL_scalar,dt,h;
+	double     vel,vx,vy,dvxdx,dvydy;
 	double     dvx[2],dvy[2];
 	double     v_gauss[2]={0.0};
@@ -533,18 +533,5 @@
 		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
 	}
-
-	//Create Artificial diffusivity once for all if requested
-	if(stabilization==1){
-		gauss=new GaussTria();
-		gauss->GaussCenter();
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
-		delete gauss;
-
-		vxaverage_input->GetInputAverage(&v_gauss[0]);
-		vyaverage_input->GetInputAverage(&v_gauss[1]);
-
-		K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
-		K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
-	}
+	h=sqrt(2*this->GetArea());
 
 	/* Start  looping on the number of gaussian points: */
@@ -591,8 +578,26 @@
 					&Ke->values[0],1);
 
-		if(stabilization==1){
+		if(stabilization==2){
+			/*Streamline upwinding*/
+			vel=sqrt(pow(vx,2.)+pow(vy,2.));
+			K[0][0]=h/(2*vel)*vx*vx;
+			K[1][0]=h/(2*vel)*vy*vx;
+			K[0][1]=h/(2*vel)*vx*vy;
+			K[1][1]=h/(2*vel)*vy*vy;
+		}
+		else if(stabilization==1){
+			/*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];
-
 			TripleMultiply( &Bprime[0][0],2,numdof,1,
 						&KDL[0][0],2,2,0,
@@ -1078,4 +1083,5 @@
 	x3=xyz_list[2][0]; y3=xyz_list[2][1];
  
+	_assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
 	return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
 }