Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 9511)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 9512)
@@ -516,5 +516,5 @@
 	int        artdiff;
 	int        i,j,ig,dim;
-	double     Jdettria ,vx,vy,dvxdx,dvydy;
+	double     Jdettria,vx,vy,dvxdx,dvydy,vel,h;
 	double     dvx[2],dvy[2];
 	double     xyz_list[NUMVERTICES][3];
@@ -546,17 +546,18 @@
 		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
 	}
-
-	/*Create Artificial diffusivity once for all if requested*/
-	if(artdiff){
-		gauss=new GaussTria();
-		gauss->GaussCenter();
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
-		delete gauss;
-
-		vxaverage_input->GetParameterAverage(&vx);
-		vyaverage_input->GetParameterAverage(&vy);
-		K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx);
-		K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy);
-	}
+	h=sqrt(2*this->GetArea());
+
+	///*Create Artificial diffusivity once for all if requested*/
+	//if(artdiff){
+	//	gauss=new GaussTria();
+	//	gauss->GaussCenter();
+	//	GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
+	//	delete gauss;
+
+	//	vxaverage_input->GetParameterAverage(&vx);
+	//	vyaverage_input->GetParameterAverage(&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:*/
@@ -596,6 +597,16 @@
 
 		if(artdiff){
+			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;
 			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,